Added new rescale method for systems

Makefile.Defines

@@ -71,7 +71,7 @@ FORTRAN = ifort
 #AR     = xiar
 
 #FORTRAN = gfortran
-#FFLAGS =  -ffree-line-length-none $(GDEBUG) $(GMEM)
+#FFLAGS =  -ffree-line-length-none $(GDEBUG) #$(GMEM)
 AR     = ar
 
 # DO NOT include in CFLAGS the "-c" option to compile object only

src/fragmentation/fragmentation.f90

@@ -314,6 +314,7 @@ subroutine calculate_system_energy(linclude_fragments)
          logical, dimension(:), allocatable :: ltmp
          logical :: lk_plpl
          class(swiftest_nbody_system), allocatable :: tmpsys
+         class(swiftest_parameters), allocatable :: tmpparam
 
          ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the
          ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by
@@ -334,27 +335,16 @@ subroutine calculate_system_energy(linclude_fragments)
                npl_new  = npl
             end if
             call setup_construct_system(tmpsys, param)
+            call tmpsys%tp%setup(0, param)
             deallocate(tmpsys%cb)
             allocate(tmpsys%cb, source=cb)
+            allocate(tmpparam, source=param)
             allocate(ltmp(npl_new))
             ltmp(:) = .false.
             ltmp(1:npl) = .true.
             call tmpsys%pl%setup(npl_new, param)
             call tmpsys%pl%fill(pl, ltmp)
-            tmpsys%pl%mass(1:npl) = tmpsys%pl%mass(1:npl) / mscale
-            tmpsys%pl%Gmass(1:npl) = tmpsys%pl%Gmass(1:npl) * tscale**2 / rscale**3
-            tmpsys%cb%mass = tmpsys%cb%mass / mscale
-            tmpsys%cb%Gmass = tmpsys%cb%Gmass * tscale**2 / rscale**3
-            tmpsys%cb%radius = tmpsys%cb%radius / rscale
-            tmpsys%pl%radius(1:npl) = tmpsys%pl%radius(1:npl) / rscale
-            tmpsys%pl%xh(:,1:npl) = tmpsys%pl%xh(:,1:npl) / rscale
-            tmpsys%pl%vh(:,1:npl) = tmpsys%pl%vh(:,1:npl) / vscale
-            tmpsys%pl%xb(:,1:npl) = tmpsys%pl%xb(:,1:npl) / rscale
-            tmpsys%pl%vb(:,1:npl) = tmpsys%pl%vb(:,1:npl) / vscale
-            tmpsys%pl%rot(:,1:npl) = tmpsys%pl%rot(:,1:npl) * tscale
-            tmpsys%cb%xb(:) = tmpsys%cb%xb(:) / rscale
-            tmpsys%cb%vb(:) = tmpsys%cb%vb(:) / vscale
-
+            call tmpsys%rescale(tmpparam, mscale, rscale, tscale)
 
             if (linclude_fragments) then ! Append the fragments if they are included
                ! Energy calculation requires the fragments to be in the system barcyentric frame, s

src/modules/swiftest_classes.f90

@@ -306,6 +306,7 @@ module swiftest_classes
       procedure :: step_spin               => tides_step_spin_system          !! Steps the spins of the massive & central bodies due to tides.
       procedure :: set_msys                => util_set_msys                   !! Sets the value of msys from the masses of system bodies.
       procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum
+      procedure :: rescale                 => util_rescale_system             !! Rescales the system into a new set of units
       procedure :: validate_ids            => util_valid_id_system            !! Validate the numerical ids passed to the system and save the maximum value
    end type swiftest_nbody_system
 
@@ -984,6 +985,13 @@ end subroutine util_fill_arr_logical
    end interface
 
    interface
+      module subroutine util_rescale_system(self, param, mscale, dscale, tscale)
+         implicit none
+         class(swiftest_nbody_system), intent(inout) :: self   !! Swiftest nbody system object
+         class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters. Returns with new values of the scale vactors and GU
+         real(DP),                     intent(in)    :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. 
+      end subroutine util_rescale_system
+
       module function util_minimize_bfgs(f, N, x0, eps, lerr) result(x1)
          use lambda_function
          implicit none
@@ -1003,38 +1011,6 @@ module subroutine util_peri_tp(self, system, param)
       end subroutine util_peri_tp
    end interface
 
-   interface util_solve_linear_system
-      module function util_solve_linear_system_d(A,b,n,lerr) result(x)
-         implicit none
-         integer(I4B),             intent(in)  :: n
-         real(DP), dimension(:,:), intent(in)  :: A
-         real(DP), dimension(:),   intent(in)  :: b
-         logical,                  intent(out) :: lerr
-         real(DP), dimension(n)                :: x
-      end function util_solve_linear_system_d
-
-      module function util_solve_linear_system_q(A,b,n,lerr) result(x)
-         implicit none
-         integer(I4B),             intent(in)  :: n
-         real(QP), dimension(:,:), intent(in)  :: A
-         real(QP), dimension(:),   intent(in)  :: b
-         logical,                  intent(out) :: lerr
-         real(QP), dimension(n)                :: x
-      end function util_solve_linear_system_q
-   end interface
-
-   interface
-      module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1)
-         use lambda_function
-         implicit none
-         class(lambda_obj),      intent(inout) :: f    !! lambda function object that has been initialized to be a function of derivatives. The object will return with components lastarg and lasteval set
-         real(DP), dimension(:), intent(in)    :: y0in !! Initial value at t=0
-         real(DP),               intent(in)    :: t1   !! Final time
-         real(DP),               intent(in)    :: dt0  !! Initial step size guess
-         real(DP),               intent(in)    :: tol  !! Tolerance on solution
-         real(DP), dimension(:), allocatable   :: y1  !! Final result
-      end function util_solve_rkf45
-   end interface
 
    interface util_resize
       module subroutine util_resize_arr_char_string(arr, nnew)
@@ -1142,6 +1118,39 @@ module subroutine util_set_rhill_approximate(self,cb)
       end subroutine util_set_rhill_approximate
    end interface
 
+   interface util_solve_linear_system
+      module function util_solve_linear_system_d(A,b,n,lerr) result(x)
+         implicit none
+         integer(I4B),             intent(in)  :: n
+         real(DP), dimension(:,:), intent(in)  :: A
+         real(DP), dimension(:),   intent(in)  :: b
+         logical,                  intent(out) :: lerr
+         real(DP), dimension(n)                :: x
+      end function util_solve_linear_system_d
+
+      module function util_solve_linear_system_q(A,b,n,lerr) result(x)
+         implicit none
+         integer(I4B),             intent(in)  :: n
+         real(QP), dimension(:,:), intent(in)  :: A
+         real(QP), dimension(:),   intent(in)  :: b
+         logical,                  intent(out) :: lerr
+         real(QP), dimension(n)                :: x
+      end function util_solve_linear_system_q
+   end interface
+
+   interface
+      module function util_solve_rkf45(f, y0in, t1, dt0, tol) result(y1)
+         use lambda_function
+         implicit none
+         class(lambda_obj),      intent(inout) :: f    !! lambda function object that has been initialized to be a function of derivatives. The object will return with components lastarg and lasteval set
+         real(DP), dimension(:), intent(in)    :: y0in !! Initial value at t=0
+         real(DP),               intent(in)    :: t1   !! Final time
+         real(DP),               intent(in)    :: dt0  !! Initial step size guess
+         real(DP),               intent(in)    :: tol  !! Tolerance on solution
+         real(DP), dimension(:), allocatable   :: y1  !! Final result
+      end function util_solve_rkf45
+   end interface
+
    interface util_sort
       module subroutine util_sort_i4b(arr)
          implicit none

src/util/util_rescale.f90

@@ -0,0 +1,51 @@
+submodule (swiftest_classes) s_util_rescale
+   use swiftest
+contains
+   module subroutine util_rescale_system(self, param, mscale, dscale, tscale)
+      !! author: David A. Minton
+      !!
+      !! Rescales an nbody system to a new set of units. Inputs are the multipliers on the mass (mscale), distance (dscale), and time units (tscale). 
+      !! Rescales all united quantities in the system, as well as the mass conversion factors, gravitational constant, and Einstein's constant in the parameter object.
+      implicit none
+      class(swiftest_nbody_system), intent(inout) :: self   !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters. Returns with new values of the scale vactors and GU
+      real(DP),                     intent(in)    :: mscale, dscale, tscale !! Scale factors for mass, distance, and time units, respectively. 
+      ! Internals
+      real(DP) :: vscale
+
+      param%MU2KG = param%MU2KG * mscale
+      param%DU2M = param%DU2M * dscale
+      param%TU2S = param%TU2S * tscale
+
+      ! Calculate the G for the system units
+      param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2))
+
+      ! Calculate the inverse speed of light in the system units
+      param%inv_c2 = einsteinC * param%TU2S / param%DU2M
+      param%inv_c2 = (param%inv_c2)**(-2)
+
+      vscale = dscale / tscale
+
+      associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody)
+
+         cb%mass = cb%mass / mscale
+         cb%Gmass = param%GU * cb%mass 
+         cb%radius = cb%radius / dscale
+         cb%xb(:) = cb%xb(:) / dscale
+         cb%vb(:) = cb%vb(:) / vscale
+         pl%mass(1:npl) = pl%mass(1:npl) / mscale
+         pl%Gmass(1:npl) = param%GU * pl%mass(1:npl) 
+         pl%radius(1:npl) = pl%radius(1:npl) / dscale
+         pl%xh(:,1:npl) = pl%xh(:,1:npl) / dscale
+         pl%vh(:,1:npl) = pl%vh(:,1:npl) / vscale
+         pl%xb(:,1:npl) = pl%xb(:,1:npl) / dscale
+         pl%vb(:,1:npl) = pl%vb(:,1:npl) / vscale
+         pl%rot(:,1:npl) = pl%rot(:,1:npl) * tscale
+
+      end associate
+
+
+      return
+   end subroutine util_rescale_system 
+
+end submodule s_util_rescale

src/util/util_resize.f90

@@ -1,6 +1,7 @@
 submodule (swiftest_classes) s_util_resize
    use swiftest
 contains
+
    module subroutine util_resize_arr_char_string(arr, nnew)
       !! author: David A. Minton
       !!