Refactored PARTICLE_FILE to PARTICLE_OUT for consistency. Added in pa…

…rticle info io to SyMBA

examples/symba_energy_momentum/param.disruption_headon.in

@@ -7,7 +7,7 @@ TP_IN          tp.in
 IN_TYPE        ASCII
 ISTEP_OUT      1                 ! output cadence every year
 BIN_OUT        bin.disruption_headon.dat
-PARTICLE_FILE  particle.disruption_headon.dat
+PARTICLE_OUT   particle.disruption_headon.dat
 OUT_TYPE       REAL8             ! double precision real output
 OUT_FORM       XV                ! osculating element output
 OUT_STAT       REPLACE

examples/symba_energy_momentum/param.disruption_off_axis.in

@@ -7,7 +7,7 @@ TP_IN          tp.in
 IN_TYPE        ASCII
 ISTEP_OUT      1                 ! output cadence every year
 BIN_OUT        bin.disruption_off_axis.dat
-PARTICLE_FILE  particle.disruption_off_axis.dat
+PARTICLE_OUT   particle.disruption_off_axis.dat
 OUT_TYPE       REAL8             ! double precision real output
 OUT_FORM       XV                ! osculating element output
 OUT_STAT       REPLACE

examples/symba_energy_momentum/param.escape.in

@@ -7,7 +7,7 @@ TP_IN          tp.in
 IN_TYPE        ASCII
 ISTEP_OUT      1                 ! output cadence every year
 BIN_OUT        bin.escape.dat
-PARTICLE_FILE  particle.escape.dat
+PARTICLE_OUT   particle.escape.dat
 OUT_TYPE       REAL8             ! double precision real output
 OUT_FORM       XV                ! osculating element output
 OUT_STAT       REPLACE

examples/symba_energy_momentum/param.sun.in

@@ -9,7 +9,7 @@ IN_TYPE        ASCII
 ISTEP_OUT      1
 ISTEP_DUMP     1
 BIN_OUT        bin.sun.dat
-PARTICLE_FILE  particle.sun.dat
+PARTICLE_OUT   particle.sun.dat
 OUT_TYPE       REAL8
 OUT_FORM       XV                ! osculating element output
 OUT_STAT       REPLACE

examples/symba_energy_momentum/param.supercatastrophic_headon.in

@@ -7,7 +7,7 @@ TP_IN          tp.in
 IN_TYPE        ASCII
 ISTEP_OUT      1                 ! output cadence every year
 BIN_OUT        bin.supercatastrophic_headon.dat
-PARTICLE_FILE  particle.supercatastrophic_headon.dat
+PARTICLE_OUT   particle.supercatastrophic_headon.dat
 OUT_TYPE       REAL8             ! double precision real output
 OUT_FORM       XV                ! osculating element output
 OUT_STAT       REPLACE

examples/symba_energy_momentum/param.supercatastrophic_off_axis.in

@@ -7,7 +7,7 @@ TP_IN          tp.in
 IN_TYPE        ASCII
 ISTEP_OUT      1                 ! output cadence every year
 BIN_OUT        bin.supercatastrophic_off_axis.dat
-PARTICLE_FILE  particle.supercatastrophic_off_axis.dat
+PARTICLE_OUT   particle.supercatastrophic_off_axis.dat
 OUT_TYPE       REAL8             ! double precision real output
 OUT_FORM       XV                ! osculating element output
 OUT_STAT       REPLACE

src/io/io.f90

@@ -498,7 +498,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
                   read(param_value, *) param%Ecollisions
                case("EUNTRACKED")
                   read(param_value, *) param%Euntracked
-               case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
+               case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_OUT", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
                case default
                   write(iomsg,*) "Unknown parameter -> ",param_name
                   iostat = -1

src/modules/symba_classes.f90

@@ -18,7 +18,7 @@ module symba_classes
    integer(I4B),          parameter :: PARTICLEUNIT     = 44 !! File unit number for the binary particle info output file
 
    type, extends(swiftest_parameters) :: symba_parameters
-      character(STRMAX)                       :: particle_file  = PARTICLE_OUTFILE !! Name of output particle information file
+      character(STRMAX)                       :: particle_out  = PARTICLE_OUTFILE !! Name of output particle information file
       real(DP)                                :: GMTINY          = -1.0_DP          !! Smallest mass that is fully gravitating
       integer(I4B), dimension(:), allocatable :: seed                              !! Random seeds
       logical                                 :: lfragmentation = .false.          !! Do fragmentation modeling instead of simple merger.
@@ -325,10 +325,11 @@ module subroutine symba_io_write_discard(self, param)
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
       end subroutine symba_io_write_discard
 
-      module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) 
+      module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) 
          implicit none
          class(symba_nbody_system),             intent(inout) :: system !! SyMBA nbody system object
-         class(symba_parameters),               intent(inout) :: param  !! Current run configuration parameters with SyMBA extensions
+         class(symba_parameters),               intent(in)    :: param  !! Current run configuration parameters with SyMBA extensions
+         logical,                     optional, intent(in)    :: lincludecb  !! Set to true to include the central body (default is false)
          integer(I4B), dimension(:),  optional, intent(in)    :: tpidx  !! Array of test particle indices to append to the particle file
          integer(I4B), dimension(:),  optional, intent(in)    :: plidx  !! Array of massive body indices to append to the particle file
       end subroutine symba_io_dump_particle_info

src/symba/symba_io.f90

@@ -2,29 +2,30 @@
    use swiftest
 contains
 
-   module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx) 
+   module subroutine symba_io_dump_particle_info(system, param, lincludecb, tpidx, plidx) 
       !! author: David A. Minton
       !!
       !! Dumps the particle information data to a file. 
       !! Pass a list of array indices for test particles (tpidx) and/or massive bodies (plidx) to append
       implicit none
       ! Arguments
       class(symba_nbody_system),             intent(inout) :: system !! SyMBA nbody system object
-      class(symba_parameters),               intent(inout) :: param  !! Current run configuration parameters with SyMBA extensions
+      class(symba_parameters),               intent(in)    :: param  !! Current run configuration parameters with SyMBA extensions
+      logical,                     optional, intent(in)    :: lincludecb  !! Set to true to include the central body (default is false)
       integer(I4B), dimension(:),  optional, intent(in)    :: tpidx  !! Array of test particle indices to append to the particle file
       integer(I4B), dimension(:),  optional, intent(in)    :: plidx  !! Array of massive body indices to append to the particle file
       ! Internals
       logical, save             :: lfirst = .true.
-      integer(I4B), parameter   :: iu = 22
+      integer(I4B), parameter   :: LUN = 22
       integer(I4B)              :: i, ierr
 
       if (.not.present(tpidx) .and. .not.present(plidx)) return
       if (lfirst) then
          select case(param%out_stat)
          case('APPEND')
-            open(unit = iu, file = param%particle_file, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr)
+            open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr)
          case('NEW', 'UNKNOWN', 'REPLACE')
-            open(unit = iu, file = param%particle_file, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr)
+            open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', iostat = ierr)
          case default
             write(*,*) 'Invalid status code',trim(adjustl(param%out_stat))
             call util_exit(FAILURE)
@@ -37,20 +38,30 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx)
 
          lfirst = .false.
       else
-         open(unit = iu, file = param%particle_file, status = 'OLD', position =  'APPEND', form = 'UNFORMATTED', iostat = ierr)
+         open(unit = LUN, file = param%particle_out, status = 'OLD', position =  'APPEND', form = 'UNFORMATTED', iostat = ierr)
          if (ierr /= 0) then
             write(*, *) "Swiftest error:"
             write(*, *) "   unable to open binary output file for APPEND"
             call util_exit(FAILURE)
          end if
       end if
 
+      if (present(lincludecb)) then
+         if (lincludecb) then
+            select type(cb => system%cb)
+            class is (symba_cb)
+               write(LUN) cb%id
+               write(LUN) cb%info
+            end select
+         end if
+      end if
+
       if (present(plidx) .and. (system%pl%nbody > 0) .and. size(plidx) > 0) then
          select type(pl => system%pl)
          class is (symba_pl)
             do i = 1, size(plidx)
-               write(iu) pl%id(plidx(i))
-               write(iu) pl%info(plidx(i))
+               write(LUN) pl%id(plidx(i))
+               write(LUN) pl%info(plidx(i))
             end do
          end select
       end if
@@ -59,13 +70,13 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx)
          select type(tp => system%tp)
          class is (symba_tp)
             do i = 1, size(tpidx)
-               write(iu) tp%id(tpidx(i))
-               write(iu) tp%info(tpidx(i))
+               write(LUN) tp%id(tpidx(i))
+               write(LUN) tp%info(tpidx(i))
             end do
          end select
       end if
 
-      close(unit = iu, iostat = ierr)
+      close(unit = LUN, iostat = ierr)
       if (ierr /= 0) then
          write(*, *) "Swiftest error:"
          write(*, *) "   unable to close particle output file"
@@ -75,6 +86,7 @@ module subroutine symba_io_dump_particle_info(system, param, tpidx, plidx)
       return
    end subroutine symba_io_dump_particle_info
 
+
    module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) 
       !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
       !!
@@ -119,6 +131,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms
                ifirst = ilast + 1
                param_value = io_get_token(line_trim, ifirst, ilast, iostat)
                select case (param_name)
+               case ("PARTICLE_OUT")
+                  param%particle_out = param_value
                case ("FRAGMENTATION")
                   call io_toupper(param_value)
                   if (param_value == "YES" .or. param_value == "T") self%lfragmentation = .true.
@@ -217,7 +231,7 @@ module subroutine symba_io_param_writer(self, unit, iotype, v_list, iostat, ioms
 
          ! Special handling is required for writing the random number seed array as its size is not known until runtime
          ! For the "SEED" parameter line, the first value will be the size of the seed array and the rest will be the seed array elements
-         write(param_name, Afmt) "PARTICLE_FILE"; write(param_value, Afmt) trim(adjustl(param%particle_file)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
+         write(param_name, Afmt) "PARTICLE_OUT"; write(param_value, Afmt) trim(adjustl(param%particle_out)); write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "GMTINY"; write(param_value, Rfmt) param%Gmtiny; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "FRAGMENTATION"; write(param_value, Lfmt)  param%lfragmentation; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          if (param%lfragmentation) then
@@ -259,7 +273,7 @@ module subroutine symba_io_read_particle(system, param)
       logical                 :: lmatch  
       type(symba_particle_info) :: tmpinfo
 
-      open(unit = LUN, file = param%particle_file, status = 'OLD', form = 'UNFORMATTED', iostat = ierr)
+      open(unit = LUN, file = param%particle_out, status = 'OLD', form = 'UNFORMATTED', iostat = ierr)
       if (ierr /= 0) then
          write(*, *) "Swiftest error:"
          write(*, *) "   unable to open binary particle file for reading"

src/symba/symba_setup.f90

@@ -12,13 +12,15 @@ module subroutine symba_setup_initialize_particle_info(system, param)
       class(symba_parameters),   intent(inout) :: param !! Current run configuration parameters with SyMBA extensions
       ! Internals
       integer(I4B) :: i
+      integer(I4B), dimension(:), allocatable :: idx
 
       select type(cb => system%cb)
       class is (symba_cb)
          cb%info%origin_type = "Central body"
          cb%info%origin_time = param%t0
          cb%info%origin_xh(:) = 0.0_DP
          cb%info%origin_vh(:) = 0.0_DP
+         call symba_io_dump_particle_info(system, param, lincludecb=.true.)
       end select
 
       select type(pl => system%pl)
@@ -29,6 +31,11 @@ module subroutine symba_setup_initialize_particle_info(system, param)
             pl%info(i)%origin_xh(:) = pl%xh(:,i)
             pl%info(i)%origin_vh(:) = pl%vh(:,i)
          end do
+         if (pl%nbody > 0) then
+            allocate(idx(pl%nbody))
+            call symba_io_dump_particle_info(system, param, plidx=[(i, i=1, pl%nbody)])
+            deallocate(idx)
+         end if
       end select
 
       select type(tp => system%tp)
@@ -38,9 +45,15 @@ module subroutine symba_setup_initialize_particle_info(system, param)
             tp%info(i)%origin_time = param%t0
             tp%info(i)%origin_xh(:) = tp%xh(:,i)
             tp%info(i)%origin_vh(:) = tp%vh(:,i)
-        end do
+         end do
+         if (tp%nbody > 0) then
+            allocate(idx(tp%nbody))
+            call symba_io_dump_particle_info(system, param, tpidx=[(i, i=1, tp%nbody)])
+            deallocate(idx)
+         end if
       end select
 
+
       return
    end subroutine symba_setup_initialize_particle_info
 

src/symba/symba_util.f90

@@ -381,7 +381,10 @@ module subroutine symba_util_rearray_pl(self, system, param)
          if (allocated(pl%xend)) deallocate(pl%xend)
 
          ! Add in any new bodies
-         call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)])
+         if (pl_adds%nbody > 0) then
+            call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)])
+            call symba_io_dump_particle_info(system, param, plidx=[(i, i = 1, pl%nbody)])
+         end if 
 
          ! If there are still bodies in the system, sort by mass in descending order and re-index
          if (pl%nbody > 0) then