Made adjustments to initial position generator in Fraggle. Also moved…

… NetCDF file I/O operations out of the nbody system constructor, as that is now used internally to make temporary systems that shouldn't be opening files.

src/fraggle/fraggle_generate.f90

@@ -252,7 +252,7 @@ module subroutine fraggle_generate_pos_vec(collider)
          ! We will treat the first two fragments of the list as special cases. 
          ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap
          loverlap(:) = .true.
-         rdistance(:) = rdistance_scale_factor * .mag.impactors%vc(:,:)
+         rdistance(:) = rdistance_scale_factor * impactors%radius(:)
          ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point
          ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be.
          call random_number(mass_rscale)
@@ -332,7 +332,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param)
       ! Internals
       real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir
       real(DP) :: v_init, v_final, mass_init, mass_final, rotmag
-      real(DP), parameter :: random_scale_factor = 0.1_DP !! The relative scale factor to apply to the random component of the rotation vector
+      real(DP), parameter :: random_scale_factor = 0.01_DP !! The relative scale factor to apply to the random component of the rotation vector
       integer(I4B) :: i
 
       associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)

src/swiftest/swiftest_util.f90

@@ -2541,13 +2541,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param)
       ! Arguments
       class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object
       class(swiftest_parameters),                intent(inout) :: param  !! Current run configuration parameters
-      ! Internals
-      type(encounter_storage)  :: encounter_history
-      type(collision_storage)  :: collision_history
 
-      allocate(swiftest_storage(param%dump_cadence) :: param%system_history)
-      allocate(swiftest_netcdf_parameters :: param%system_history%nc)
-      call param%system_history%reset()
 
       select case(param%integrator)
       case (INT_BS)
@@ -2602,32 +2596,6 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param)
             allocate(symba_list_plpl :: nbody_system%plpl_encounter)
             allocate(collision_list_plpl :: nbody_system%plpl_collision)
 
-            if (param%lenc_save_trajectory .or. param%lenc_save_closest) then
-               allocate(encounter_netcdf_parameters :: encounter_history%nc)
-               call encounter_history%reset()
-               select type(nc => encounter_history%nc)
-               class is (encounter_netcdf_parameters)
-                  nc%file_name = ENCOUNTER_OUTFILE
-                  if (.not.param%lrestart) then
-                     call nc%initialize(param)
-                     call nc%close()
-                  end if
-               end select
-               allocate(nbody_system%encounter_history, source=encounter_history)
-            end if
-
-            allocate(collision_netcdf_parameters :: collision_history%nc)
-            call collision_history%reset()
-            select type(nc => collision_history%nc)
-            class is (collision_netcdf_parameters)
-               nc%file_name = COLLISION_OUTFILE
-               if (.not.param%lrestart) then
-                  call nc%initialize(param)
-                  call nc%close()
-               end if
-            end select
-            allocate(nbody_system%collision_history, source=collision_history)
-
          end select
       case (INT_RINGMOONS)
          write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled'
@@ -2638,15 +2606,6 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param)
 
       allocate(swiftest_particle_info :: nbody_system%cb%info)
 
-      select case(param%collision_model)
-      case("MERGE")
-         allocate(collision_basic :: nbody_system%collider)
-      case("BOUNCE")
-         allocate(collision_bounce :: nbody_system%collider)
-      case("FRAGGLE")
-         allocate(collision_fraggle :: nbody_system%collider)
-      end select
-      call nbody_system%collider%setup(nbody_system)
 
       nbody_system%t = param%tstart
 
@@ -2698,8 +2657,11 @@ module subroutine swiftest_util_setup_initialize_system(self, param)
       class(swiftest_nbody_system), intent(inout) :: self   !! Swiftest nbody_system object
       class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters
 
-      associate(nbody_system => self, cb => self%cb, pl => self%pl, tp => self%tp, nc => param%system_history%nc)
+      allocate(swiftest_storage(param%dump_cadence) :: param%system_history)
+      allocate(swiftest_netcdf_parameters :: param%system_history%nc)
+      call param%system_history%reset()
 
+      associate(nbody_system => self, cb => self%cb, pl => self%pl, tp => self%tp, nc => param%system_history%nc)
          call nbody_system%read_in(param)
          call nbody_system%validate_ids(param)
          call nbody_system%set_msys()
@@ -2722,7 +2684,6 @@ module subroutine swiftest_util_setup_initialize_system(self, param)
          nc%file_name = param%outfile
          call nbody_system%write_frame(param) 
          call nc%close()
-
       end associate
 
       return

src/symba/symba_module.f90

@@ -22,7 +22,6 @@ module symba
    real(DP),     private, parameter :: RHSCALE = 6.5_DP
    real(DP),     private, parameter :: RSHELL  = 0.48075_DP
 
-
    !> SyMBA central body particle class
    type, extends(helio_cb) :: symba_cb
    end type symba_cb

src/symba/symba_util.f90

@@ -9,6 +9,7 @@
 
 submodule(symba) s_symba_util
    use swiftest
+   use fraggle
 contains
 
    module subroutine symba_util_append_pl(self, source, lsource_mask)
@@ -267,13 +268,53 @@ module subroutine symba_util_setup_initialize_system(self, param)
       ! Arguments
       class(symba_nbody_system),  intent(inout) :: self    !! SyMBA nbody_system object
       class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
+      ! Internals
+      type(encounter_storage)  :: encounter_history
+      type(collision_storage)  :: collision_history
 
       ! Call parent method
       associate(nbody_system => self)
          call helio_util_setup_initialize_system(nbody_system, param)
          call nbody_system%pltp_encounter%setup(0_I8B)
          call nbody_system%plpl_encounter%setup(0_I8B)
          call nbody_system%plpl_collision%setup(0_I8B)
+
+         if (param%lenc_save_trajectory .or. param%lenc_save_closest) then
+            allocate(encounter_netcdf_parameters :: encounter_history%nc)
+            call encounter_history%reset()
+            select type(nc => encounter_history%nc)
+            class is (encounter_netcdf_parameters)
+               nc%file_name = ENCOUNTER_OUTFILE
+               if (.not.param%lrestart) then
+                  call nc%initialize(param)
+                  call nc%close()
+               end if
+            end select
+            allocate(nbody_system%encounter_history, source=encounter_history)
+         end if
+
+         allocate(collision_netcdf_parameters :: collision_history%nc)
+         call collision_history%reset()
+         select type(nc => collision_history%nc)
+         class is (collision_netcdf_parameters)
+            nc%file_name = COLLISION_OUTFILE
+            if (.not.param%lrestart) then
+               call nc%initialize(param)
+               call nc%close()
+            end if
+         end select
+         allocate(nbody_system%collision_history, source=collision_history)
+
+         select case(param%collision_model)
+         case("MERGE")
+            allocate(collision_basic :: nbody_system%collider)
+         case("BOUNCE")
+            allocate(collision_bounce :: nbody_system%collider)
+         case("FRAGGLE")
+            allocate(collision_fraggle :: nbody_system%collider)
+         end select
+         call nbody_system%collider%setup(nbody_system)
+
       end associate
 
       return