! Copyright 2025 - David Minton
! This file is part of Swiftest.
! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License 
! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
! You should have received a copy of the GNU General Public License along with Swiftest. 
! If not, see: https://www.gnu.org/licenses. 

module swiftest
   !! author: David A. Minton
   !!
   !! This module serves to combine all of the Swiftest project modules under a single umbrella so that they can be accessed from 
   !! individual submodule implementations with a simple "use swiftest" line.
   !! 
   !! The project structure is divided into a heirarchy of modules. The lowest level of the heirarchy are the modules called in the
   !! "use" statements below. Next the "swiftest" modules (this one), and finally each individual integrator (and potential future 
   !! integrators) sit at the top. This structure is a consequence of two competing constraints:
   !! 1) The desire that much of the basic functionality of the code is modular, such that new functionality can be easily added 
   !!    without altering too much of the basic code.
   !! 2) Adhering to Modern Fortran's typing rules.
   !!  
   !! A set of "base" types is defined in the base module. These define classes of objects, (i.e. central body, massive body, and 
   !! test particles) and other major types used throughout the project. However, none of the derived data types are defined with 
   !! concrete type-bound procedures attached (only abstract procedures). However, the *interfaces* of type-bound procedures are 
   !! defined using the base types as arguments. Because of the typing rules of Modern Fortran's type-bound procedure overrides, any
   !! non-pass arguments(i.e. arguments not named self) must be identical in all extended types. Because some of the basic 
   !! functionality in the project is split across multiple modules, we cannot define type-bound procedures in base class objects 
   !! until the all interfaces are defined. In order to avoid these dependency issues and not end up with a massive base class with
   !! every possibly type-bound procedure interface in the project (thus reducing the modularity of the project), the type-bound 
   !! procedures are added to the base types here. 
   !!
   !! Structuring this code this way adds somewhat to the verbosity of the code. The main thing that has to happen is that for any 
   !! procedures where one wishes to make use of an type-bound procedures defined for arguments at the swiftest-type level or 
   !! higher, but that are passsed to base-level procedures, must have their arguments wrapped in a select type(...); class is(...) 
   !! construct in order to "reveal" the procedures. This is done throughout the project at the beginning of many procedures (along 
   !! with copious amounts of associate(...) statements, in order to help with code readibility)
   !!
   !! Adapted from David E. Kaufmann's Swifter routine: module_swifter.f90
   use globals
   use operators
   use lambda_function
   use base
   use encounter
   use collision
   use walltime
   use io_progress_bar
   use netcdf_io
   use solver
   !use advisor_annotate
   !$ use omp_lib
   implicit none
   public

   type, extends(netcdf_parameters) :: swiftest_netcdf_parameters
   contains
      procedure :: initialize => swiftest_io_netcdf_initialize_output 
         !! Initialize a set of parameters used to identify a NetCDF output object
      procedure :: get_valid_masks => swiftest_io_netcdf_get_valid_masks   
         !! Gets logical masks indicating which bodies are valid pl and tp type at the current time
      procedure :: open            => swiftest_io_netcdf_open 
         !! Opens a NetCDF file and does the variable inquiries to activate variable ids
      procedure :: flush           => swiftest_io_netcdf_flush 
         !! Flushes a NetCDF file by closing it then opening it again
   end type swiftest_netcdf_parameters


   type, extends(base_storage) :: swiftest_storage
      class(swiftest_netcdf_parameters), allocatable :: nc 
         !! NetCDF object attached to this storage object
   contains
      procedure :: dump             => swiftest_io_dump_storage 
         !! Dumps storage object contents to file
      procedure :: dealloc          => swiftest_util_dealloc_storage
         !! Resets a storage object by deallocating all items and resetting the frame counter to 0
      procedure :: get_index_values => swiftest_util_get_vals_storage
         !! Gets the unique values of the indices of a storage object (i.e. body id or time value)
      procedure :: make_index_map   => swiftest_util_index_map_storage
         !! Maps body id values to storage index values so we don't have to use unlimited dimensions for id
      procedure :: take_snapshot    => swiftest_util_snapshot_system 
         !! Takes a snapshot of the nbody_system for later file storage
      final     ::                     swiftest_final_storage
   end type swiftest_storage


   ! The following extended types or their children should be used, where possible, as the base of any types defined in additional 
   ! modules, such as new integrators. 
   type, extends(base_parameters) :: swiftest_parameters
   contains
      procedure :: dump        => swiftest_io_dump_param
      procedure :: reader      => swiftest_io_param_reader
      procedure :: writer      => swiftest_io_param_writer
      procedure :: read_in     => swiftest_io_read_in_param
      procedure :: set_display => swiftest_io_set_display_param
   end type swiftest_parameters


   !> Class definition for the kinship relationships used in bookkeeping multiple collisions bodies in a single time step.
   type, extends(base_kinship) :: swiftest_kinship
      integer(I4B)                            :: parent 
         !! Index of parent particle
      integer(I4B)                            :: nchild
         !! number of children in merger list
      integer(I4B), dimension(:), allocatable :: child
         !! Index of children particles
   contains
      procedure :: dealloc  => swiftest_util_dealloc_kin 
         !! Deallocates all allocatable arrays
      final :: swiftest_final_kin
         !! Finalizes the Swiftest kinship object - deallocates all allocatables
   end type swiftest_kinship


   type, extends(base_particle_info) :: swiftest_particle_info
      character(len=NAMELEN)    :: name            
         !! Non-unique name
      character(len=NAMELEN)    :: particle_type   
         !! String containing a description of the particle type (e.g. Central Body,  Massive Body, Test Particle)
      character(len=NAMELEN)    :: origin_type     
         !! String containing a description of the origin of the particle (e.g. Initial  Conditions, Disruption, etc.)
      real(DP)                  :: origin_time  = -huge(1.0_DP)
         !! The time of the particle's formation
      integer(I4B)              :: collision_id = 0   
         !! The ID of the collision that formed the particle
      real(DP), dimension(NDIM) :: origin_rh  = [0.0_DP, 0.0_DP, 0.0_DP] 
         !! The heliocentric distance vector at the time of the particle's formation
      real(DP), dimension(NDIM) :: origin_vh  = [0.0_DP, 0.0_DP, 0.0_DP] 
         !! The heliocentric velocity vector at the time of the particle's formation
      real(DP)                  :: discard_time = huge(1.0_DP) 
         !! The time of the particle's discard
      character(len=NAMELEN)    :: status          
         !! Particle status description: Active, Merged, Fragmented, etc.
      real(DP), dimension(NDIM) :: discard_rh = [0.0_DP, 0.0_DP, 0.0_DP]     
         !! The heliocentric distance vector at the time of the particle's discard
      real(DP), dimension(NDIM) :: discard_vh = [0.0_DP, 0.0_DP, 0.0_DP] 
         !! The heliocentric velocity vector at the time of the particle's discard
      integer(I4B)              :: discard_body_id = -1
         !! The id of the other body involved in the discard (0 if no other body involved)
   contains
      procedure :: copy      => swiftest_util_copy_particle_info  
         !! Copies one set of information object components into another,  component-by-component
      procedure :: set_value => swiftest_util_set_particle_info   
         !! Sets one or more values of the particle information metadata object
   end type swiftest_particle_info


   !> An abstract class for a generic collection of Swiftest bodies
   type, abstract, extends(base_object) :: swiftest_body
         !! Superclass that defines the generic elements of a Swiftest particle 
      logical                                                    :: lfirst = .true. 
         !! Run the current step as a first
      integer(I4B)                                               :: nbody = 0       
         !! Number of bodies
      integer(I4B),                  dimension(:),   allocatable :: id              
         !! Identifier 
      type(swiftest_particle_info),  dimension(:),   allocatable :: info            
         !! Particle metadata information
      logical,                       dimension(:),   allocatable :: lmask           
         !! Logical mask used to select a subset of bodies when performing certain operations (drift, kick, accel, etc.)
      integer(I4B),                  dimension(:),   allocatable :: status          
         !! An integrator-specific status indicator 
      logical,                       dimension(:),   allocatable :: ldiscard        
         !! Body should be discarded
      logical,                       dimension(:),   allocatable :: lcollision      
         !! flag indicating whether body has merged with another this time step
      logical,                       dimension(:),   allocatable :: lencounter      
         !! flag indicating whether body is part of an encounter this time step
      real(DP),                      dimension(:),   allocatable :: mu              
         !! G * (Mcb + [m])
      real(DP),                      dimension(:,:), allocatable :: rh              
         !! Heliocentric position
      real(DP),                      dimension(:,:), allocatable :: vh              
         !! Heliocentric velocity
      real(DP),                      dimension(:,:), allocatable :: rb              
         !! Barycentric position
      real(DP),                      dimension(:,:), allocatable :: vb              
         !! Barycentric velocity
      real(DP),                      dimension(:,:), allocatable :: ah              
         !! Total heliocentric acceleration
      real(DP),                      dimension(:,:), allocatable :: aobl            
         !! Barycentric accelerations of bodies due to central body oblatenes
      real(DP),                      dimension(:,:), allocatable :: agr             
         !! Acceleration due to post-Newtonian correction
      real(DP),                      dimension(:,:), allocatable :: atide           
         !! Tanngential component of acceleration of bodies due to tides
      real(DP),                      dimension(:),   allocatable :: ir3h            
         !! Inverse heliocentric radius term (1/rh**3)
      integer(I4B),                  dimension(:),   allocatable :: isperi          
         !! perihelion passage flag
      real(DP),                      dimension(:),   allocatable :: peri            
         !! perihelion distance
      real(DP),                      dimension(:),   allocatable :: atp             
         !! semimajor axis following perihelion passage
      real(DP),                      dimension(:),   allocatable :: a               
         !! Semimajor axis (pericentric distance for a parabolic orbit)
      real(DP),                      dimension(:),   allocatable :: e               
         !! Eccentricity
      real(DP),                      dimension(:),   allocatable :: inc             
         !! Inclination
      real(DP),                      dimension(:),   allocatable :: capom           
         !! Longitude of ascending node
      real(DP),                      dimension(:),   allocatable :: omega           
         !! Argument of pericenter
      real(DP),                      dimension(:),   allocatable :: capm            
         !! Mean anomaly
      ! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the
      ! component list, such as setup_body and util_spill
   contains
      procedure(abstract_discard_body), deferred :: discard
      procedure(abstract_kick_body),    deferred :: kick     
      procedure(abstract_set_mu),       deferred :: set_mu
      procedure(abstract_step_body),    deferred :: step
      procedure(abstract_accel),        deferred :: accel

      ! These are concrete because the implementation is the same for all types of particles
      procedure :: drift           => swiftest_drift_body                   
         !! Loop through bodies and call Danby drift routine on heliocentric variables
      procedure :: v2pv            => swiftest_gr_vh2pv_body                
         !! Converts from velocity to psudeovelocity for GR calculations using symplectic integrators
      procedure :: pv2v            => swiftest_gr_pv2vh_body                
         !! Converts from psudeovelocity to velocity for GR calculations using symplectic integrators
      procedure :: read_frame_bin  => swiftest_io_read_frame_body           
         !! I/O routine for writing out a single frame of time-series data for the central body
      procedure :: read_in         => swiftest_io_read_in_body              
         !! Read in body initial conditions from an ascii file
      procedure :: write_frame     => swiftest_io_netcdf_write_frame_body   
         !! I/O routine for writing out a single frame of time-series data for all bodies in the nbody_system in NetCDF format  
      procedure :: write_info      => swiftest_io_netcdf_write_info_body    
         !! Dump contents of particle information metadata to file
      procedure :: el2xv           => swiftest_orbel_el2xv_vec              
         !! Convert orbital elements to position and velocity vectors
      procedure :: xv2el           => swiftest_orbel_xv2el_vec              
         !! Convert position and velocity vectors to orbital elements 
      procedure :: setup           => swiftest_util_setup_body              
         !! A constructor that sets the number of bodies and allocates all allocatable arrays
      procedure :: accel_user      => swiftest_user_kick_getacch_body       
         !! Add user-supplied heliocentric accelerations to planets
      procedure :: append          => swiftest_util_append_body             
         !! Appends elements from one structure to another
      procedure :: dealloc         => swiftest_util_dealloc_body            
         !! Deallocates all allocatable arrays
      procedure :: fill            => swiftest_util_fill_body               
         !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
      procedure :: get_peri        => swiftest_util_peri_body               
         !! Determine nbody_system pericenter passages for test particles 
      procedure :: resize          => swiftest_util_resize_body             
         !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small.
      procedure :: set_ir3         => swiftest_util_set_ir3h                
         !! Sets the inverse heliocentric radius term (1/rh**3)
      procedure :: sort            => swiftest_util_sort_body               
         !! Sorts body arrays by a sortable componen
      procedure :: rearrange       => swiftest_util_sort_rearrange_body     
         !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods
      procedure :: spill           => swiftest_util_spill_body              
         !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
      procedure :: save_discard => swiftest_util_save_discard_body
         !! Saves a snapshot of the this body to the collision storage object
      generic   :: read_frame      => read_frame_bin                        
         !! Add the generic read frame for Fortran binary files
   end type swiftest_body


   type, abstract, extends(base_object) :: swiftest_cb
   !> An abstract class for a generic central body in a Swiftest simulation
      class(swiftest_particle_info), allocatable  :: info              
         !! Particle metadata information
      integer(I4B)                                :: id       = 0      
         !! External identifier (unique)
      real(DP)                                    :: mass     = 0.0_DP 
         !! Central body mass (units MU)
      real(DP)                                    :: Gmass    = 0.0_DP 
         !! Central mass gravitational term G * mass (units GU * MU)
      real(DP)                                    :: radius   = 0.0_DP 
         !! Central body radius (units DU)
      real(DP)                                    :: density  = 1.0_DP 
         !! Central body mass density - calculated internally (units MU / DU**3)
      real(DP)                                    :: j2rp2    = 0.0_DP 
         !! J2*R^2 term for central body
      real(DP)                                    :: j4rp4    = 0.0_DP 
         !! J4*R^4 term for central body
      real(DP), dimension(:,:,:), allocatable     :: c_lm              
         !! Spherical Harmonics coefficients for the central body
      real(DP), dimension(NDIM)                   :: aobl     = 0.0_DP 
         !! Barycentric acceleration due to central body oblatenes
      real(DP), dimension(NDIM)                   :: atide    = 0.0_DP 
         !! Barycentric acceleration due to central body oblatenes
      real(DP), dimension(NDIM)                   :: aoblbeg  = 0.0_DP 
         !! Barycentric acceleration due to central body oblatenes at beginning of step
      real(DP), dimension(NDIM)                   :: aoblend  = 0.0_DP 
         !! Barycentric acceleration due to central body oblatenes at end of step
      real(DP), dimension(NDIM)                   :: atidebeg = 0.0_DP 
         !! Barycentric acceleration due to central body oblatenes at beginning of step
      real(DP), dimension(NDIM)                   :: atideend = 0.0_DP 
         !! Barycentric acceleration due to central body oblatenes at end of step
      real(DP), dimension(NDIM)                   :: rb       = 0.0_DP 
         !! Barycentric position (units DU)
      real(DP), dimension(NDIM)                   :: vb       = 0.0_DP 
         !! Barycentric velocity (units DU / TU)
      real(DP), dimension(NDIM)                   :: agr      = 0.0_DP 
         !! Acceleration due to post-Newtonian correction
      real(DP), dimension(NDIM)                   :: Ip       = 0.0_DP 
         !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. 
      real(DP), dimension(NDIM)                   :: rot      = 0.0_DP 
         !! Body rotation vector in inertial coordinate frame (units rad / TU)
      real(DP)                                    :: rotphase = 0.0_DP 
         !! Body rotation phase about the rotation pole (0 to 1)
      real(DP)                                    :: k2       = 0.0_DP 
         !! Tidal Love number
      real(DP)                                    :: Q        = 0.0_DP 
         !! Tidal quality factor
      real(DP)                                    :: tlag     = 0.0_DP 
         !! Tidal phase lag angle
      real(DP), dimension(NDIM)                   :: L0       = 0.0_DP 
         !! Initial angular momentum of the central body
      real(DP), dimension(NDIM)                   :: dL       = 0.0_DP 
         !! Change in angular momentum of the central body
      real(DP)                                    :: GM0      = 0.0_DP 
         !! Initial G*mass of the central body
      real(DP)                                    :: dGM      = 0.0_DP 
         !! Change in G*mass of the central body
      real(DP)                                    :: R0       = 0.0_DP 
         !! Initial radius of the central body
      real(DP)                                    :: dR       = 0.0_DP 
         !! Change in the radius of the central body
   contains
      procedure :: rotphase_update => swiftest_drift_cb_rotphase_update 
         !! updates the central body rotation phase
      procedure :: dealloc      => swiftest_util_dealloc_cb          
         !! Deallocates all allocatables and resets all values to defaults 
      procedure :: read_in      => swiftest_io_read_in_cb            
         !! Read in central body initial conditions from an ASCII file
      procedure :: write_frame  => swiftest_io_netcdf_write_frame_cb 
         !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format  
      procedure :: write_info   => swiftest_io_netcdf_write_info_cb  
         !! Dump contents of particle information metadata to file
   end type swiftest_cb


   type, abstract, extends(swiftest_body) :: swiftest_pl
         !! Superclass that defines the generic elements of a Swiftest particle 
      real(DP),                dimension(:),   allocatable :: mass    
         !! Body mass (units MU)
      real(DP),                dimension(:),   allocatable :: Gmass   
         !! Mass gravitational term G * mass (units GU * MU)
      real(DP),                dimension(:),   allocatable :: rhill   
         !! Body mass (units MU)
      real(DP),                dimension(:),   allocatable :: renc    
         !! Critical radius for close encounters
      real(DP),                dimension(:),   allocatable :: radius  
         !! Body radius (units DU)
      real(DP),                dimension(:),   allocatable :: density 
         !! Body mass density - calculated internally (units MU / DU**3)
      real(DP),                dimension(:,:), allocatable :: rbeg    
         !! Position at beginning of step
      real(DP),                dimension(:,:), allocatable :: rend    
         !! Position at end of step
      real(DP),                dimension(:,:), allocatable :: vbeg    
         !! Velocity at beginning of step
      real(DP),                dimension(:,:), allocatable :: Ip      
         !! Unitless principal moments of inertia (I1, I2, I3) / (MR**2). Principal axis rotation assumed. 
      real(DP),                dimension(:,:), allocatable :: rot     
         !! Body rotation vector in inertial coordinate frame (units rad / TU)
      real(DP),                dimension(:),   allocatable :: k2      
         !! Tidal Love number
      real(DP),                dimension(:),   allocatable :: Q       
         !! Tidal quality factor
      real(DP),                dimension(:),   allocatable :: tlag    
         !! Tidal phase lag
      integer(I4B),            dimension(:,:), allocatable :: k_plpl  
         !! Index array used to convert flattened the body-body comparison upper triangular matrix
      integer(I8B)                                         :: nplpl   
         !! Number of body-body comparisons in the flattened upper triangular matrix
      type(swiftest_kinship),  dimension(:),   allocatable :: kin        
         !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step
      logical,                 dimension(:),   allocatable :: lmtiny     
         !! flag indicating whether this body is below the GMTINY cutoff value
      integer(I4B)                                         :: nplm = 0   
         !! number of bodies above the GMTINY limit
      integer(I8B)                                         :: nplplm     
         !! Number of body (all massive)-body (only those above GMTINY) comparisons in the flattened upper triangular matrix 
      integer(I4B),            dimension(:),   allocatable :: nplenc     
         !! number of encounters with other planets this time step
      integer(I4B),            dimension(:),   allocatable :: ntpenc     
         !! number of encounters with test particles this time step

      ! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the
      !    component list, such as setup_pl and util_spill_pl
   contains
      ! Massive body-specific concrete methods 
      ! These are concrete because they are the same implemenation for all integrators
      procedure :: make_impactors => swiftest_util_make_impactors_pl 
         !! Make impactors out of the current kinship relationships
      procedure :: discard => swiftest_discard_pl             
         !! Placeholder method for discarding massive bodies 
      procedure :: accel_int => swiftest_kick_getacch_int_pl    
         !! Compute direct cross (third) term heliocentric accelerations of massive bodies
      procedure :: accel_non_spherical_cb => swiftest_non_spherical_cb_acc_pl             
         !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
      procedure :: setup => swiftest_util_setup_pl          
         !! A base constructor that sets the number of bodies and allocates and initializes all arrays  
      ! procedure :: accel_tides    => tides_kick_getacch_pl           
         !! Compute the accelerations of bodies due to tidal interactions with the central body
      procedure :: append => swiftest_util_append_pl         
         !! Appends elements from one structure to another
      procedure :: h2b  => swiftest_util_coord_h2b_pl      
         !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity)
      procedure :: b2h => swiftest_util_coord_b2h_pl      
         !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity)
      procedure :: vh2vb  => swiftest_util_coord_vh2vb_pl    
         !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
      procedure :: vb2vh          => swiftest_util_coord_vb2vh_pl    
         !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
      procedure :: rh2rb          => swiftest_util_coord_rh2rb_pl    
         !! Convert massive bodies from heliocentric to barycentric coordinates (position only)
      procedure :: dealloc        => swiftest_util_dealloc_pl        
         !! Deallocates all allocatable arrays
      procedure :: fill           => swiftest_util_fill_pl           
         !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
      procedure :: flatten        => swiftest_util_flatten_eucl_plpl 
         !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
      procedure :: rearray        => swiftest_util_rearray_pl        
         !! Clean up the massive body structures to remove discarded bodies and add new bodies
      procedure :: resize         => swiftest_util_resize_pl         
         !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small.
      procedure :: reset_kinship  => swiftest_util_reset_kinship_pl  
         !! Resets the kinship status of bodies
      procedure :: set_beg_end    => swiftest_util_set_beg_end_pl    
         !! Sets the beginning and ending positions and velocities of planets.
      procedure :: set_mu         => swiftest_util_set_mu_pl         
         !! Method used to construct the vectorized form of the central body mass
      procedure :: set_rhill      => swiftest_util_set_rhill         
         !! Calculates the Hill's radii for each body
      procedure :: set_renc_I4B   => swiftest_util_set_renc_I4B      
         !! Sets the critical radius for encounter given an inpput integer scale factor
      procedure :: set_renc_DP    => swiftest_util_set_renc_DP       
         !! Sets the critical radius for encounter given an input real scale factor
      procedure :: sort           => swiftest_util_sort_pl           
         !! Sorts body arrays by a sortable component
      procedure :: rearrange      => swiftest_util_sort_rearrange_pl 
         !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods
      procedure :: spill          => swiftest_util_spill_pl          
         !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
      generic   :: set_renc       => set_renc_I4B, set_renc_DP 
   end type swiftest_pl


   type, abstract, extends(swiftest_body) :: swiftest_tp
      
         !! Superclass that defines the generic elements of a Swiftest test particle 
      integer(I4B), dimension(:,:), allocatable :: k_pltp 
         !! Index array used to convert flattened the body-body comparison upper triangular matrix
      integer(I8B)                              :: npltp  
         !! Number of pl-tp comparisons in the flattened upper triangular matrix
      integer(I4B), dimension(:),   allocatable :: nplenc 
         !! number of encounters with planets this time step
      
      ! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the
      ! component list, such as swiftest_util_setup_tp and util_spill_tp
   contains
      ! Test particle-specific concrete methods 
      ! These are concrete because they are the same implemenation for all integrators
      procedure :: discard   => swiftest_discard_tp             
         !! Check to see if test particles should be discarded based on their positions relative to the massive bodies
      procedure :: accel_int => swiftest_kick_getacch_int_tp    
         !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies
      procedure :: accel_non_spherical_cb => swiftest_non_spherical_cb_acc_tp             
         !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
      procedure :: setup     => swiftest_util_setup_tp          
         !! A base constructor that sets the number of bodies and 
      procedure :: append    => swiftest_util_append_tp         
         !! Appends elements from one structure to another
      procedure :: h2b       => swiftest_util_coord_h2b_tp      
         !! Convert test particles from heliocentric to barycentric coordinates (position and velocity)
      procedure :: b2h       => swiftest_util_coord_b2h_tp      
         !! Convert test particles from barycentric to heliocentric coordinates (position and velocity)
      procedure :: vb2vh     => swiftest_util_coord_vb2vh_tp    
         !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
      procedure :: vh2vb     => swiftest_util_coord_vh2vb_tp    
         !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
      procedure :: rh2rb     => swiftest_util_coord_rh2rb_tp    
         !! Convert test particles from heliocentric to barycentric coordinates (position only)
      procedure :: dealloc   => swiftest_util_dealloc_tp        
         !! Deallocates all allocatable arrays
      procedure :: fill      => swiftest_util_fill_tp           
         !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic)
      procedure :: rearray   => swiftest_util_rearray_tp        
         !! Clean up the test particle structures to remove discarded bodies 
      procedure :: resize    => swiftest_util_resize_tp         
         !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small.
      procedure :: set_mu    => swiftest_util_set_mu_tp         
         !! Method used to construct the vectorized form of the central body mass
      procedure :: sort      => swiftest_util_sort_tp           
         !! Sorts body arrays by a sortable component
      procedure :: rearrange => swiftest_util_sort_rearrange_tp 
         !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods
      procedure :: spill     => swiftest_util_spill_tp          
         !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
   end type swiftest_tp


   !> An abstract class for a basic Swiftest nbody system 
   type, abstract, extends(base_nbody_system) :: swiftest_nbody_system
         !! This superclass contains a minimial nbody_system of a set of test particles (tp), massive bodies (pl), and a central 
         !! body (cb). 
      class(swiftest_cb),         allocatable :: cb                
         !! Central body data structure
      class(swiftest_pl),         allocatable :: pl                
         !! Massive body data structure
      class(swiftest_tp),         allocatable :: tp                
         !! Test particle data structure
      class(swiftest_tp),         allocatable :: tp_discards       
         !! Discarded test particle data structure
      class(swiftest_pl),         allocatable :: pl_discards       
         !! Discarded massive body particle data structure
      class(swiftest_pl),         allocatable :: pl_adds           
         !! List of added bodies in mergers or collisions
      class(swiftest_tp),         allocatable :: tp_adds           
         !! List of added bodies in mergers or collisions
      class(encounter_list),      allocatable :: pltp_encounter    
         !! List of massive body-test particle encounters in a single step
      class(encounter_list),      allocatable :: plpl_encounter    
         !! List of massive body-massive body encounters in a single step
      class(collision_list_plpl), allocatable :: plpl_collision    
         !! List of massive body-massive body collisions in a single step
      class(collision_list_pltp), allocatable :: pltp_collision    
         !! List of massive body-test particle collisions in a single step
      class(collision_basic),     allocatable :: collider          
         !! Collision system object
      class(encounter_storage),   allocatable :: encounter_history 
         !! Stores encounter history for later retrieval and saving to file
      class(collision_storage),   allocatable :: collision_history 
         !! Stores encounter history for later retrieval and saving to file
      integer(I4B)                    :: maxid = -1             
         !! The current maximum particle id number 
      real(DP)                        :: t = -1.0_DP            
         !! Integration current time
      real(DP)                        :: GMtot = 0.0_DP         
         !! Total nbody_system mass - used for barycentric coordinate conversion
      real(DP)                        :: ke_orbit = 0.0_DP      
         !! nbody_system orbital kinetic energy
      real(DP)                        :: ke_rot = 0.0_DP       
         !! nbody_system rotational kinetic energy
      real(DP)                        :: pe = 0.0_DP            
         !! nbody_system potential energy
      real(DP)                        :: be = 0.0_DP            
         !! nbody_system binding energy of all bodies
      real(DP)                        :: te = 0.0_DP            
         !! nbody_system total energy
      real(DP)                        :: oblpot = 0.0_DP        
         !! nbody_system potential energy due to oblateness of the central body
      real(DP), dimension(NDIM)       :: L_orbit = 0.0_DP       
         !! nbody_system orbital angular momentum vector
      real(DP), dimension(NDIM)       :: L_rot = 0.0_DP        
         !! nbody_system rotational angular momentum vector
      real(DP), dimension(NDIM)       :: L_total = 0.0_DP       
         !! nbody_system angular momentum vector
      real(DP)                        :: ke_orbit_orig = 0.0_DP 
         !! Initial orbital kinetic energy
      real(DP)                        :: ke_rot_orig = 0.0_DP  
         !! Initial rotational kinetic energy
      real(DP)                        :: pe_orig = 0.0_DP       
         !! Initial potential energy
      real(DP)                        :: be_orig = 0.0_DP       
         !! Initial gravitational binding energy
      real(DP)                        :: te_orig = 0.0_DP       
         !! Initial total energy (sum of all sources of energy tracked)
      real(DP)                        :: be_cb   = 0.0_DP       
         !! Binding energy of central body (usually orders of magnitude larger than the rest of the system, and therefore tracked 
         !! seperately)
      real(DP)                        :: E_orbit_orig = 0.0_DP  
         !! Initial orbital energy
      real(DP)                        :: GMtot_orig = 0.0_DP    
         !! Initial nbody_system mass
      real(DP), dimension(NDIM)       :: L_total_orig = 0.0_DP  
         !! Initial total angular momentum vector
      real(DP), dimension(NDIM)       :: L_orbit_orig = 0.0_DP  
         !! Initial orbital angular momentum
      real(DP), dimension(NDIM)       :: L_rot_orig = 0.0_DP   
         !! Initial rotational angular momentum vector
      real(DP), dimension(NDIM)       :: L_escape = 0.0_DP      
         !! Angular momentum of bodies that escaped the nbody_system (used for bookeeping)
      real(DP)                        :: GMescape = 0.0_DP      
         !! Mass of bodies that escaped the nbody_system (used for bookeeping)
      real(DP)                        :: E_collisions = 0.0_DP  
         !! Energy lost from nbody_system due to collisions
      real(DP)                        :: E_untracked = 0.0_DP   
         !! Energy gained from nbody_system due to escaped bodies

      ! Energy, momentum, and mass errors (used in error reporting)
      real(DP)                        :: ke_orbit_error    = 0.0_DP
      real(DP)                        :: ke_rot_error     = 0.0_DP
      real(DP)                        :: pe_error          = 0.0_DP
      real(DP)                        :: be_error          = 0.0_DP
      real(DP)                        :: E_orbit_error     = 0.0_DP
      real(DP)                        :: Ecoll_error       = 0.0_DP
      real(DP)                        :: E_untracked_error = 0.0_DP
      real(DP)                        :: te_error          = 0.0_DP
      real(DP)                        :: L_orbit_error     = 0.0_DP
      real(DP)                        :: L_rot_error      = 0.0_DP
      real(DP)                        :: L_escape_error    = 0.0_DP
      real(DP)                        :: L_total_error     = 0.0_DP
      real(DP)                        :: Mtot_error        = 0.0_DP
      real(DP)                        :: Mescape_error     = 0.0_DP

      logical                         :: lbeg                 
         !! True if this is the beginning of a step. This is used so that test particle steps can be calculated separately from 
         !! massive bodies.  Massive body variables are saved at half steps, and passed to the test particles
      logical                         :: lfirst_io   = .true.  
         !! Flag to indicate that this is the first time to write to a file
      logical                         :: lfirst_peri = .true.  
         !! Flag to indicate that this is the first pericenter passage
   contains
      !> Each integrator will have its own version of the step
      procedure(abstract_step_system), deferred :: step

      ! Concrete classes that are common to the basic integrator (only test particles considered for discard)
      procedure :: discard                 => swiftest_discard_system                              
         !! Perform a discard step on the nbody_system
      procedure :: compact_output          => swiftest_io_compact_output                           
         !! Prints out out terminal output when display_style is set to COMPACT
      procedure :: conservation_report     => swiftest_io_conservation_report                      
         !! Compute energy and momentum and print out the change with time
      procedure :: display_run_information => swiftest_io_display_run_information                  
         !! Displays helpful information about the run
      procedure :: dump                    => swiftest_io_dump_system                              
         !! Dump the state of the nbody_system to a file
      procedure :: get_t0_values           => swiftest_io_netcdf_get_t0_values_system              
         !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output.
      procedure :: read_frame              => swiftest_io_netcdf_read_frame_system                 
         !! Read in a frame of input data from file
      procedure :: read_hdr                => swiftest_io_netcdf_read_hdr_system                   
         !! Read a header for an output frame in NetCDF format
      procedure :: write_hdr               => swiftest_io_netcdf_write_hdr_system                  
         !! Write a header for an output frame in NetCDF format
      procedure :: read_particle_info      => swiftest_io_netcdf_read_particle_info_system         
         !! Read in particle metadata from file
      procedure :: read_in                 => swiftest_io_read_in_system                           
         !! Reads the initial conditions for an nbody system
      procedure :: write_frame             => swiftest_io_netcdf_write_frame_system                
         !! Write a frame of input data from file
      procedure :: obl_pot                 => swiftest_obl_pot_system                              
         !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body
      procedure :: dealloc                 => swiftest_util_dealloc_system                         
         !! Deallocates all allocatables and resets all values to defaults. Acts as a base for a finalizer
      procedure :: get_energy_and_momentum => swiftest_util_get_energy_and_momentum_system         
         !! Calculates the total nbody_system energy and momentum
      procedure :: get_idvals              => swiftest_util_get_idvalues_system                    
         !! Returns an array of all id values in use in the nbody_system
      procedure :: rescale                 => swiftest_util_rescale_system                         
         !! Rescales the nbody_system into a new set of units
      procedure :: initialize_output_file  => swiftest_io_initialize_output_file_system                  
         !! Write a frame of input data from file
      procedure :: initialize              => swiftest_util_setup_initialize_system                
         !! Initialize the nbody_system from input files
      procedure :: init_particle_info      => swiftest_util_setup_initialize_particle_info_system  
         !! Initialize the nbody_system from input files
    ! procedure :: step_rot               => tides_step_rot_system                               
         !! Steps the rotations of the massive & central bodies due to tides.
      procedure :: set_msys                => swiftest_util_set_msys                               
         !! Sets the value of msys from the masses of nbody_system bodies.
      procedure :: validate_ids            => swiftest_util_valid_id_system                        
         !! Validate the numerical ids passed to the nbody_system and save the maximum value
#ifdef COARRAY
      procedure :: coarray_collect         => swiftest_coarray_collect_system                      
         !! Collects all the test particles from other images into the image #1 test particle system
      procedure :: coarray_distribute      => swiftest_coarray_distribute_system                   
         !! Distributes test particles from image #1 out to all images.
      procedure :: coarray_balance         => swiftest_coarray_balance_system                      
         !! Checks whether or not the test particle coarrays need to be rebalanced.
#endif
   end type swiftest_nbody_system


   abstract interface

      subroutine abstract_accel(self, nbody_system, param, t, lbeg)
         import swiftest_body, swiftest_nbody_system, swiftest_parameters, DP
         class(swiftest_body),         intent(inout) :: self   
            !! Swiftest body data structure
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),   intent(inout) :: param  
            !! Current run configuration parameters 
         real(DP),                     intent(in)    :: t      
            !! Current simulation time
         logical,                      intent(in)    :: lbeg   
            !! Optional argument that determines whether or not this is the beginning or end of the step
      end subroutine abstract_accel

      subroutine abstract_discard_body(self, nbody_system, param) 
         import swiftest_body, swiftest_nbody_system, swiftest_parameters
         class(swiftest_body),              intent(inout) :: self   
            !! Swiftest body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(inout) :: param  
            !! Current run configuration parameters 
      end subroutine abstract_discard_body

      subroutine abstract_kick_body(self, nbody_system, param, t, dt, lbeg)
         import swiftest_body, swiftest_nbody_system, swiftest_parameters, DP
         implicit none
         class(swiftest_body),              intent(inout) :: self   
            !! Swiftest generic body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system objec
         class(swiftest_parameters),        intent(inout) :: param  
            !! Current run configuration parameters 
         real(DP),                          intent(in)    :: t      
            !! Current time
         real(DP),                          intent(in)    :: dt     
            !! Stepsize
         logical,                           intent(in)    :: lbeg   
            !! Logical flag indicating whether this is the beginning of the half step or not. 
      end subroutine abstract_kick_body

      subroutine abstract_set_mu(self, cb) 
         import swiftest_body, swiftest_cb
         class(swiftest_body), intent(inout) :: self 
            !! Swiftest body object
         class(swiftest_cb),   intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine abstract_set_mu

      subroutine abstract_step_body(self, nbody_system, param, t, dt)
         import DP, swiftest_body, swiftest_nbody_system, swiftest_parameters
         implicit none
         class(swiftest_body),         intent(inout) :: self         
            !! Swiftest body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody_system object
         class(swiftest_parameters),   intent(inout) :: param        
            !! Current run configuration parameters 
         real(DP),                     intent(in)    :: t            
            !! Simulation time
         real(DP),                     intent(in)    :: dt           
            !! Current stepsize
      end subroutine abstract_step_body

      subroutine abstract_step_system(self, param, t, dt)
         import DP, swiftest_nbody_system, swiftest_parameters
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self  
            !! Swiftest nbody_system object
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
         real(DP),                          intent(in)    :: t     
            !! Simulation time
         real(DP),                          intent(in)    :: dt    
            !! Current stepsize
      end subroutine abstract_step_system
   end interface

   interface
      module subroutine swiftest_discard_pl(self, nbody_system, param)
         implicit none
         class(swiftest_pl),           intent(inout) :: self   
            !! Swiftest massive body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),   intent(inout) :: param  
            !! Current run configuration parameter
      end subroutine swiftest_discard_pl

      module subroutine swiftest_discard_system(self, param)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self  
            !! Swiftest nbody_system object
         class(swiftest_parameters),   intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_discard_system

      module subroutine swiftest_discard_tp(self, nbody_system, param)
         implicit none
         class(swiftest_tp),           intent(inout) :: self   
            !! Swiftest test particle object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),   intent(inout) :: param  
            !! Current run configuration parameters
      end subroutine swiftest_discard_tp

      module subroutine swiftest_drift_all(mu, x, v, n, param, dt, lmask, iflag)
         implicit none
         real(DP), dimension(:),     intent(in)    :: mu    
            !! Vector of gravitational constants
         real(DP), dimension(:,:),   intent(inout) :: x, v  
            !! Position and velocity vectors
         integer(I4B),               intent(in)    :: n     
            !! number of bodies
         class(swiftest_parameters),     intent(in)    :: param 
            !! Current run configuration parameters
         real(DP),                   intent(in)    :: dt    
            !! Stepsize
         logical, dimension(:),      intent(in)    :: lmask 
            !! Logical mask of size self%nbody that determines which bodies to drift.
         integer(I4B), dimension(:), intent(out)   :: iflag 
            !! Vector of error flags. 0 means no problem
      end subroutine swiftest_drift_all

      module subroutine swiftest_drift_body(self, nbody_system, param, dt)
         implicit none
         class(swiftest_body),              intent(inout) :: self   
            !! Swiftest particle data structure
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(in)    :: param  
            !! Current run configuration parameters
         real(DP),                          intent(in)    :: dt     
            !! Stepsize
      end subroutine swiftest_drift_body

      pure elemental module subroutine swiftest_drift_one(mu, rx, ry, rz, vx, vy, vz, dt, iflag)
         !$omp declare simd(swiftest_drift_one)
         implicit none
         real(DP),     intent(in)       :: mu    
            !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift
         real(DP),     intent(inout)    :: rx, ry, rz, vx, vy, vz  
            !! Position and velocity of body to drift
         real(DP),     intent(in)       :: dt    
            !! Step size
         integer(I4B), intent(out)      :: iflag 
            !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
      end subroutine swiftest_drift_one

      module subroutine swiftest_drift_cb_rotphase_update(self, param, dt)
         implicit none
         class(swiftest_cb),           intent(inout) :: self   
            !! Swiftest central body data structure
         class(swiftest_parameters),   intent(in)    :: param  
            !! Current run configuration parameters 
         real(DP),                     intent(in)    :: dt     
            !! Stepsize
      end subroutine

      module subroutine swiftest_driver(integrator, param_file_name, display_style)
         implicit none
         character(len=:), intent(in), allocatable :: integrator      
            !! Symbolic code of the requested integrator  
         character(len=:), intent(in), allocatable :: param_file_name 
            !! Name of the input parameters file
         character(len=:), intent(in), allocatable :: display_style   
            !! {"PROGRESS", "CLASSIC", "QUIET", "COMPACT"}). 
            !! Default is "PROGRESS"   
      end subroutine swiftest_driver

      pure module subroutine swiftest_gr_kick_getaccb_ns_body(self, nbody_system, param)
         implicit none
         class(swiftest_body),              intent(inout) :: self   
            !! Swiftest generic body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(in)    :: param  
            !! Current run configuration parameters 
      end subroutine swiftest_gr_kick_getaccb_ns_body

      pure module subroutine swiftest_gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) 
         implicit none
         real(DP), dimension(:),     intent(in)  :: mu     
            !! Gravitational constant
         real(DP), dimension(:,:),   intent(in)  :: x      
            !! Position vectors
         logical,  dimension(:),     intent(in)  :: lmask  
            !! Logical mask indicating which bodies to compute
         integer(I4B),               intent(in)  :: n      
            !! Total number of bodies
         real(DP),                   intent(in)  :: inv_c2 
            !! Inverse speed of light squared: 1 / c**2
         real(DP), dimension(:,:),   intent(out) :: agr    
            !! Accelerations
      end subroutine swiftest_gr_kick_getacch

      pure elemental module subroutine swiftest_gr_p4_pos_kick(inv_c2, rx, ry, rz, vx, vy, vz, dt)
         implicit none
         real(DP),  intent(in)    :: inv_c2     
            !! One over speed of light squared (1/c**2)
         real(DP),  intent(inout) :: rx, ry, rz 
            !! Position vector
         real(DP),  intent(in)    :: vx, vy, vz 
            !! Velocity vector
         real(DP),  intent(in)    :: dt         
            !! Step size
      end subroutine swiftest_gr_p4_pos_kick

      pure module subroutine swiftest_gr_pseudovel2vel(param, mu, rh, pv, vh) 
         implicit none
         class(swiftest_parameters), intent(in)  :: param 
            !! Current run configuration parameters 
         real(DP),                   intent(in)  :: mu    
            !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
         real(DP), dimension(:),     intent(in)  :: rh    
            !! Swiftestcentric position vector 
         real(DP), dimension(:),     intent(in)  :: pv    
            !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
         real(DP), dimension(:),     intent(out) :: vh    
            !! Swiftestcentric velocity vector 
      end subroutine swiftest_gr_pseudovel2vel

      pure module subroutine swiftest_gr_pv2vh_body(self, param)
         implicit none
         class(swiftest_body),       intent(inout) :: self  
            !! Swiftest particle object
         class(swiftest_parameters), intent(in)    :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_gr_pv2vh_body

      pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv)
         implicit none
         class(swiftest_parameters), intent(in)  :: param 
            !! Current run configuration parameters 
         real(DP),                   intent(in)  :: mu    
            !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
         real(DP), dimension(:),     intent(in)  :: rh    
            !! Swiftestcentric position vector 
         real(DP), dimension(:),     intent(in)  :: vh    
            !! Swiftestcentric velocity vector 
         real(DP), dimension(:),     intent(out) :: pv    
            !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
      end subroutine swiftest_gr_vel2pseudovel

      pure module subroutine swiftest_gr_vh2pv_body(self, param)
         implicit none
         class(swiftest_body),       intent(inout) :: self  
            !! Swiftest particle object
         class(swiftest_parameters), intent(in)    :: param 
            !! Current run configuration parameters
      end subroutine swiftest_gr_vh2pv_body

      module subroutine swiftest_io_compact_output(self, param, timer)
         implicit none
         class(swiftest_nbody_system), intent(in) :: self  
            !! Swiftest nbody system object   
         class(swiftest_parameters),        intent(in) :: param 
            !! Input colleciton of user-defined parameters
         class(*),                          intent(in) :: timer 
            !! Object used for computing elapsed wall time
      end subroutine swiftest_io_compact_output

      module subroutine swiftest_io_conservation_report(self, param, lterminal)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self      
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(inout) :: param     
            !! Input colleciton of user-defined parameters
         logical,                           intent(in)    :: lterminal 
            !! Indicates whether to output information to the terminal screen
      end subroutine swiftest_io_conservation_report

      module subroutine swiftest_io_display_run_information(self, param, integration_timer, phase)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self 
            !! Swiftest nbody system object
         class(swiftest_parameters),   intent(inout) :: param 
            !! Current run configuration parameters 
         type(walltimer),              intent(inout) :: integration_timer 
            !! Object used for computing elapsed wall time
         character(len=*), optional,   intent(in)    :: phase 
            !! One of "first" or "last" 
      end subroutine swiftest_io_display_run_information

      module subroutine swiftest_io_dump_param(self, param_file_name)
         implicit none
         class(swiftest_parameters),intent(in)    :: self            
            !! Output collection of parameters
         character(len=*),          intent(in)    :: param_file_name 
            !! Parameter input file name (i.e. param.in)
      end subroutine swiftest_io_dump_param

      module subroutine swiftest_io_dump_system(self, param, system_history)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self   
            !! Swiftest nbody_system object
         class(swiftest_parameters),   intent(inout) :: param  
            !! Current run configuration parameters 
         class(swiftest_storage),      intent(inout) :: system_history    
            !! Stores the system history between output dumps
      end subroutine swiftest_io_dump_system

      module subroutine swiftest_io_dump_storage(self, param)
         implicit none
         class(swiftest_storage), intent(inout) :: self   
            !! Swiftest simulation history storage object
         class(swiftest_parameters), intent(inout) :: param  
            !! Current run configuration parameters 
      end subroutine swiftest_io_dump_storage

      module subroutine swiftest_io_get_args(integrator, param_file_name, display_style, from_cli) 
         implicit none
         character(len=:), allocatable, intent(inout) :: integrator      
            !! Symbolic code of the requested integrator  
         character(len=:), allocatable, intent(inout) :: param_file_name 
            !! Name of the input parameters file
         character(len=:), allocatable, intent(inout) :: display_style   
            !! {"PROGRESS", "CLASSIC", "QUIET", "COMPACT"}. 
            !! Default is "PROGRESS"   
         logical,                       intent(in)    :: from_cli        
            !! If true, get command-line arguments. Otherwise, use the values of the input variables
      end subroutine swiftest_io_get_args

      module function swiftest_io_get_token(buffer, ifirst, ilast, ierr) result(token)
         implicit none
         character(len=*), intent(in)    :: buffer         
            !! Input string buffer
         integer(I4B),     intent(inout) :: ifirst         
            !! Index of the buffer at which to start the search for a token
         integer(I4B),     intent(out)   :: ilast          
            !! Index of the buffer at the end of the returned token
         integer(I4B),     intent(out)   :: ierr           
            !! Error code
         character(len=:), allocatable   :: token          
            !! Returned token string
      end function swiftest_io_get_token

      module subroutine swiftest_io_log_one_message(file, message)
         implicit none
         character(len=*), intent(in) :: file   
            !! Name of file to log
         character(len=*), intent(in) :: message
      end subroutine swiftest_io_log_one_message
   
      module subroutine swiftest_io_log_start(param, file, header)
         implicit none
         class(swiftest_parameters), intent(in) :: param  
            !! Current Swiftest run configuration parameters
         character(len=*),           intent(in) :: file   
            !! Name of file to log
         character(len=*),           intent(in) :: header 
            !! Header to print at top of log file
      end subroutine swiftest_io_log_start

      module subroutine swiftest_io_netcdf_flush(self, param)
         implicit none
         class(swiftest_netcdf_parameters), intent(inout) :: self  
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_netcdf_flush

      module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param) 
         implicit none
         class(swiftest_nbody_system),      intent(inout) :: self  
            !! Swiftest nbody system object
         class(swiftest_netcdf_parameters), intent(inout) :: nc     
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_netcdf_get_t0_values_system

      module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask, plmmask, Gmtiny)
         implicit none
         class(swiftest_netcdf_parameters),  intent(inout)          :: self    
            !! Parameters used to identify a particular NetCDF dataset
         logical, dimension(:), allocatable, intent(out)            :: plmask  
            !! Logical mask indicating which bodies are massive bodies
         logical, dimension(:), allocatable, intent(out)            :: tpmask  
            !! Logical mask indicating which bodies are test particles
         logical, dimension(:), allocatable, intent(out), optional  :: plmmask 
            !! Logical mask indicating which bodies are fully interacting massive bodies
         real(DP),                           intent(in),  optional  :: Gmtiny  
            !! The cutoff G*mass between semi-interacting and fully interacting massive bodies
      end subroutine swiftest_io_netcdf_get_valid_masks

      module subroutine swiftest_io_netcdf_initialize_output(self, param)
         implicit none
         class(swiftest_netcdf_parameters), intent(inout) :: self  
            !! Parameters used to for writing a NetCDF dataset to file
         class(swiftest_parameters),        intent(in)    :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_netcdf_initialize_output

      module subroutine swiftest_io_netcdf_open(self, param, readonly)
         implicit none
         class(swiftest_netcdf_parameters), intent(inout) :: self     
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param    
            !! Current run configuration parameters
         logical, optional,                 intent(in)    :: readonly 
            !! Logical flag indicating that this should be open read only
      end subroutine swiftest_io_netcdf_open

      module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ierr)
         implicit none
         class(swiftest_nbody_system),      intent(inout) :: self  
            !! Swiftest nbody_system object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to for reading a NetCDF dataset to file
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
         integer(I4B)                                     :: ierr  
            !! Error code: returns 0 if the read is successful
      end function swiftest_io_netcdf_read_frame_system

      module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) 
         implicit none
         class(swiftest_nbody_system),      intent(inout) :: self  
            !! Swiftest nbody system object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to for reading a NetCDF dataset to file
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine swiftest_io_netcdf_read_hdr_system

      module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, plmask, tpmask)
         implicit none
         class(swiftest_nbody_system),      intent(inout) :: self   
            !! Swiftest nbody system object
         class(swiftest_netcdf_parameters), intent(inout) :: nc     
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param  
            !! Current run configuration parameters
         logical, dimension(:),             intent(in)    :: plmask 
            !! Logical array indicating which index values belong to massive bodies
         logical, dimension(:),             intent(in)    :: tpmask 
            !! Logical array indicating which index values belong to test particles
      end subroutine swiftest_io_netcdf_read_particle_info_system

      module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param)
         implicit none
         class(swiftest_body),              intent(in)    :: self  
            !! Swiftest base object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to for writing a NetCDF dataset to file
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_netcdf_write_frame_body

      module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param)
         implicit none
         class(swiftest_cb),                intent(inout)    :: self  
            !! Swiftest base object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to for writing a NetCDF dataset to file
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_netcdf_write_frame_cb

      module subroutine swiftest_io_netcdf_write_frame_system(self, nc, param)
         implicit none
         class(swiftest_nbody_system),      intent(inout) :: self  
            !! Swiftest nbody_system object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to for writing a NetCDF dataset to file
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_netcdf_write_frame_system

      module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) 
         implicit none
         class(swiftest_nbody_system),      intent(in)    :: self  
            !! Swiftest nbody system object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to for writing a NetCDF dataset to file
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine swiftest_io_netcdf_write_hdr_system

      module subroutine swiftest_io_netcdf_write_info_body(self, nc, param)
         implicit none
         class(swiftest_body),              intent(in)    :: self  
            !! Swiftest particle object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine swiftest_io_netcdf_write_info_body

      module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param)
         implicit none
         class(swiftest_cb),                intent(in)    :: self  
            !! Swiftest particle object
         class(swiftest_netcdf_parameters), intent(inout) :: nc    
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine swiftest_io_netcdf_write_info_cb

      module subroutine swiftest_io_remove_nul_char(string)
         implicit none
         character(len=*), intent(inout) :: string 
            !! String to remove nul characters from
      end subroutine swiftest_io_remove_nul_char

      module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) 
         implicit none
         class(swiftest_parameters), intent(inout) :: self       
            !! Collection of parameters
         integer(I4B),               intent(in)    :: unit       
            !! File unit number
         character(len=*),           intent(in)    :: iotype     
            !! Dummy argument passed to the  input/output procedure contains the text from the char-literal-constant, prefixed 
            !! with DT. If you do not include a char-literal-constant, the iotype argument contains only DT.
         character(len=*),           intent(in)    :: v_list(:)  
            !! The first element passes the integrator code to the reader
         integer(I4B),               intent(out)   :: iostat     
            !! IO status code
         character(len=*),           intent(inout) :: iomsg      
            !! Message to pass if iostat /= 0
      end subroutine swiftest_io_param_reader

      module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, iomsg) 
         implicit none
         class(swiftest_parameters), intent(in)    :: self      
            !! Collection of parameters
         integer(I4B),               intent(in)    :: unit      
            !! File unit number
         character(len=*),           intent(in)    :: iotype    
            !! Dummy argument passed to the  input/output procedure contains the text from the char-literal-constant, prefixed 
            !! with DT. If you do not include a char-literal-constant, the iotype argument contains only DT.
         integer(I4B),               intent(in)    :: v_list(:) 
            !! Not used in this procedure
         integer(I4B),               intent(out)   :: iostat    
            !! IO status code
         character(len=*),           intent(inout) :: iomsg     
            !! Message to pass if iostat /= 0
      end subroutine swiftest_io_param_writer
   end interface

   interface io_param_writer_one
      module subroutine swiftest_io_param_writer_one_char(param_name, param_value, unit)
         implicit none
         character(len=*), intent(in)    :: param_name  
            !! Name of parameter to print
         character(len=*), intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),     intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_char

      module subroutine swiftest_io_param_writer_one_DP(param_name, param_value, unit)
         implicit none
         character(len=*), intent(in)    :: param_name  
            !! Name of parameter to print
         real(DP),         intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),     intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_DP

      module subroutine swiftest_io_param_writer_one_DParr(param_name, param_value, unit)
         implicit none
         character(len=*),       intent(in)    :: param_name  
            !! Name of parameter to print
         real(DP), dimension(:), intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),           intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_DParr

      module subroutine swiftest_io_param_writer_one_I4B(param_name, param_value, unit)
         implicit none
         character(len=*), intent(in)    :: param_name  
            !! Name of parameter to print
         integer(I4B),     intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),     intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_I4B

      module subroutine swiftest_io_param_writer_one_I4Barr(param_name, param_value, unit)
         implicit none
         character(len=*),           intent(in)    :: param_name  
            !! Name of parameter to print
         integer(I4B), dimension(:), intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),               intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_I4Barr

      module subroutine swiftest_io_param_writer_one_I8B(param_name, param_value, unit)
         implicit none
         character(len=*), intent(in)    :: param_name  
            !! Name of parameter to print
         integer(I8B),     intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),     intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_I8B

      module subroutine swiftest_io_param_writer_one_logical(param_name, param_value, unit)
         implicit none
         character(len=*), intent(in)    :: param_name  
            !! Name of parameter to print
         logical,          intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),     intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_logical

#ifdef QUADPREC
      module subroutine swiftest_io_param_writer_one_QP(param_name, param_value, unit)
         implicit none
         character(len=*), intent(in)    :: param_name  
            !! Name of parameter to print
         real(QP),         intent(in)    :: param_value 
            !! Value of parameter to print
         integer(I4B),     intent(in)    :: unit        
            !! Open file unit number to print parameter to
      end subroutine swiftest_io_param_writer_one_QP
#endif
   end interface io_param_writer_one

   interface

      module subroutine swiftest_io_read_in_body(self,param)
         implicit none
         class(swiftest_body),       intent(inout) :: self  
            !! Swiftest base object
         class(swiftest_parameters), intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_read_in_body

      module subroutine swiftest_io_read_in_cb(self,param)
         implicit none
         class(swiftest_cb),         intent(inout) :: self  
            !! Swiftest base object
         class(swiftest_parameters), intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_read_in_cb

      module subroutine swiftest_io_read_in_param(self, param_file_name) 
         implicit none
         class(swiftest_parameters), intent(inout) :: self            
            !! Current run configuration parameters
         character(len=*),           intent(in)    :: param_file_name 
            !! Parameter input file name (i.e. param.in)
      end subroutine swiftest_io_read_in_param

      module subroutine swiftest_io_read_in_system(self, nc, param)
         implicit none
         class(swiftest_nbody_system),      intent(inout) :: self
         class(swiftest_netcdf_parameters), intent(inout) :: nc     
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param
      end subroutine swiftest_io_read_in_system

      module function swiftest_io_read_frame_body(self, iu, param) result(ierr)
         implicit none
         class(swiftest_body),       intent(inout) :: self  
            !! Swiftest body object
         integer(I4B),               intent(inout) :: iu    
            !! Unit number for the output file to write frame to
         class(swiftest_parameters), intent(inout) :: param 
            !! Current run configuration parameters 
         integer(I4B)                              :: ierr  
            !! Error code: returns 0 if the read is successful
      end function swiftest_io_read_frame_body

      module function swiftest_io_read_frame_system(self, iu, param) result(ierr)
         implicit none
         class(swiftest_nbody_system),intent(inout) :: self  
            !! Swiftest nbody_system object
         integer(I4B),                     intent(inout) :: iu    
            !! Unit number for the output file to read frame from
         class(swiftest_parameters),       intent(inout) :: param 
            !! Current run configuration parameters 
         integer(I4B)                               :: ierr  
            !! Error code: returns 0 if the read is successful
      end function swiftest_io_read_frame_system

      module subroutine swiftest_io_set_display_param(self, display_style)
         implicit none
         class(swiftest_parameters), intent(inout) :: self            
            !! Current run configuration parameters
         character(*),               intent(in)    :: display_style   
            !! Style of the output display 
      end subroutine swiftest_io_set_display_param

      module subroutine swiftest_io_toupper(string)
         implicit none
         character(*), intent(inout) :: string 
            !! String to make upper case
      end subroutine swiftest_io_toupper

      module subroutine swiftest_io_initialize_output_file_system(self, nc, param)
         implicit none
         class(swiftest_nbody_system),      intent(inout) :: self   
            !! Swiftest nbody_system object
         class(swiftest_netcdf_parameters), intent(inout) :: nc     
            !! Parameters used to identify a particular NetCDF dataset
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters 
      end subroutine swiftest_io_initialize_output_file_system

      module subroutine swiftest_kick_getacch_int_pl(self, param)
         implicit none
         class(swiftest_pl),         intent(inout) :: self  
            !! Swiftest massive body object
         class(swiftest_parameters), intent(inout) :: param 
            !! Current swiftest run configuration parameters
      end subroutine swiftest_kick_getacch_int_pl

      module subroutine swiftest_kick_getacch_int_tp(self, param, GMpl, rhp, npl)
         implicit none
         class(swiftest_tp),         intent(inout) :: self  
            !! Swiftest test particle object
         class(swiftest_parameters), intent(inout) :: param 
            !! Current swiftest run configuration parameters
         real(DP), dimension(:),     intent(in)    :: GMpl  
            !! Massive body masses
         real(DP), dimension(:,:),   intent(in)    :: rhp   
            !! Massive body position vectors
         integer(I4B),               intent(in)    :: npl   
            !! Number of active massive bodies
      end subroutine swiftest_kick_getacch_int_tp
   end interface

   interface swiftest_kick_getacch_int_all
      module subroutine swiftest_kick_getacch_int_all_flat_rad_pl(npl, nplpl, k_plpl, r, Gmass, radius, acc)
         implicit none
         integer(I4B),                 intent(in)             :: npl    
            !! Number of massive bodies
         integer(I8B),                 intent(in)             :: nplpl  
            !! Number of massive body interactions to compute
         integer(I4B), dimension(:,:), intent(in)             :: k_plpl 
            !! Array of interaction pair indices (flattened upper triangular matrix)
         real(DP),     dimension(:,:), intent(in)             :: r      
            !! Position vector array
         real(DP),     dimension(:),   intent(in)             :: Gmass  
            !! Array of massive body G*mass
         real(DP),     dimension(:),   intent(in)             :: radius 
            !! Array of massive body radii
         real(DP),     dimension(:,:), intent(inout)          :: acc    
            !! Acceleration vector array 
      end subroutine swiftest_kick_getacch_int_all_flat_rad_pl

      module subroutine swiftest_kick_getacch_int_all_flat_norad_pl(npl, nplpl, k_plpl, r, Gmass, acc)
         implicit none
         integer(I4B),                 intent(in)             :: npl    
            !! Number of massive bodies
         integer(I8B),                 intent(in)             :: nplpl  
            !! Number of massive body interactions to compute
         integer(I4B), dimension(:,:), intent(in)             :: k_plpl 
            !! Array of interaction pair indices (flattened upper triangular matrix)
         real(DP),     dimension(:,:), intent(in)             :: r      
            !! Position vector array
         real(DP),     dimension(:),   intent(in)             :: Gmass  
            !! Array of massive body G*mass
         real(DP),     dimension(:,:), intent(inout)          :: acc    
            !! Acceleration vector array 
      end subroutine swiftest_kick_getacch_int_all_flat_norad_pl

      module subroutine swiftest_kick_getacch_int_all_tri_rad_pl(npl, nplm, r, Gmass, radius, acc)
         implicit none
         integer(I4B),                 intent(in)             :: npl    
            !! Total number of massive bodies
         integer(I4B),                 intent(in)             :: nplm   
            !! Number of fully interacting massive bodies
         real(DP),     dimension(:,:), intent(in)             :: r      
            !! Position vector array
         real(DP),     dimension(:),   intent(in)             :: Gmass  
            !! Array of massive body G*mass
         real(DP),     dimension(:),   intent(in)             :: radius 
            !! Array of massive body radii
         real(DP),     dimension(:,:), intent(inout)          :: acc    
            !! Acceleration vector array 
      end subroutine swiftest_kick_getacch_int_all_tri_rad_pl

      module subroutine swiftest_kick_getacch_int_all_tri_norad_pl(npl, nplm, r, Gmass, acc)
         implicit none
         integer(I4B),                 intent(in)             :: npl    
            !! Total number of massive bodies
         integer(I4B),                 intent(in)             :: nplm   
            !! Number of fully interacting massive bodies
         real(DP),     dimension(:,:), intent(in)             :: r      
            !! Position vector array
         real(DP),     dimension(:),   intent(in)             :: Gmass  
            !! Array of massive body G*mass
         real(DP),     dimension(:,:), intent(inout)          :: acc    
            !! Acceleration vector array 
      end subroutine swiftest_kick_getacch_int_all_tri_norad_pl

      module subroutine swiftest_kick_getacch_int_all_tp(ntp, npl, rtp, rpl, GMpl, lmask, acc)
         implicit none
         integer(I4B),                 intent(in)    :: ntp   
            !! Number of test particles
         integer(I4B),                 intent(in)    :: npl   
            !! Number of massive bodies
         real(DP),     dimension(:,:), intent(in)    :: rtp   
            !! Test particle position vector array
         real(DP),     dimension(:,:), intent(in)    :: rpl   
            !! Massive body particle position vector array
         real(DP),     dimension(:),   intent(in)    :: GMpl  
            !! Array of massive body G*mass
         logical,      dimension(:),   intent(in)    :: lmask 
            !! Logical mask indicating which test particles should be computed
         real(DP),     dimension(:,:), intent(inout) :: acc   
            !! Acceleration vector array 
      end subroutine swiftest_kick_getacch_int_all_tp
   end interface

   interface
      pure module subroutine swiftest_kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
         !$omp declare simd(swiftest_kick_getacch_int_one_pl)
         implicit none
         real(DP), intent(in)  :: rji2            
            !! Square of distance between the two bodies
         real(DP), intent(in)  :: xr, yr, zr      
            !! Distances between the two bodies in x, y, and z directions
         real(DP), intent(in)  :: Gmi             
            !! G*mass of body i
         real(DP), intent(in)  :: Gmj             
            !! G*mass of body j
         real(DP), intent(inout) :: axi, ayi, azi 
            !! Acceleration vector components of body i
         real(DP), intent(inout) :: axj, ayj, azj 
            !! Acceleration vector components of body j
      end subroutine swiftest_kick_getacch_int_one_pl

      pure module subroutine swiftest_kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, az)
         !$omp declare simd(swiftest_kick_getacch_int_one_tp)
         implicit none
         real(DP), intent(in)  :: rji2         
            !! Square of distance between the test particle and massive body
         real(DP), intent(in)  :: xr, yr, zr   
            !! Distances between the two bodies in x, y, and z directions
         real(DP), intent(in)  :: Gmpl         
            !! G*mass of massive body
         real(DP), intent(inout) :: ax, ay, az 
            !! Acceleration vector components of test particle
      end subroutine swiftest_kick_getacch_int_one_tp

      module subroutine swiftest_obl_rot_matrix(n, rot, rot_matrix, rot_matrix_inv)
         implicit none
         integer(I4B),             intent(in)           :: n              
            !! Number of bodies
         real(DP), dimension(NDIM), intent(in)       :: rot               
            !! Central body rotation vector
         real(DP), dimension(NDIM, NDIM), intent(inout) :: rot_matrix     
            !! rotation matrix 
         real(DP), dimension(NDIM, NDIM), intent(inout) :: rot_matrix_inv 
            !! inverse of the rotation matrix
      end subroutine swiftest_obl_rot_matrix

      module subroutine swiftest_obl_acc(n, GMcb, j2rp2, j4rp4, rh, lmask, aobl, rot, GMpl, aoblcb)
         implicit none
         integer(I4B),             intent(in)            :: n      
            !! Number of bodies
         real(DP),                 intent(in)            :: GMcb   
            !! Central body G*Mass
         real(DP),                 intent(in)            :: j2rp2  
            !! J2 * R**2 for the central body
         real(DP),                 intent(in)            :: j4rp4  
            !! J4 * R**4 for the central body
         real(DP), dimension(:,:), intent(in)            :: rh     
            !! Heliocentric positions of bodies
         logical,  dimension(:),   intent(in)            :: lmask  
            !! Logical mask of bodies to compute aobl
         real(DP), dimension(:,:), intent(out)           :: aobl   
            !! Barycentric acceleration of bodies due to central body oblateness
         real(DP), dimension(NDIM), intent(in)        :: rot    
            !! Central body rotation matrix
         real(DP), dimension(:),   intent(in),  optional :: GMpl   
            !! Masses of input bodies if they are not test particles
         real(DP), dimension(:),   intent(out), optional :: aoblcb 
            !! Barycentric acceleration of central body (only needed if input bodies are massive)
      end subroutine swiftest_obl_acc

      module subroutine swiftest_non_spherical_cb_acc_pl(self, nbody_system)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
      end subroutine swiftest_non_spherical_cb_acc_pl

      module subroutine swiftest_non_spherical_cb_acc_tp(self, nbody_system)
         implicit none
         class(swiftest_tp), intent(inout) :: self  
            !! Swiftest test particle object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
      end subroutine swiftest_non_spherical_cb_acc_tp

      module subroutine swiftest_obl_pot_system(self)
         implicit none
         class(swiftest_nbody_system), intent(inout)  :: self   
            !! Swiftest nbody system object
      end subroutine swiftest_obl_pot_system

      pure elemental module subroutine swiftest_orbel_el2xv(mu, a, ie, inc_deg, capom_deg, omega_deg, capm_deg, &
                                                                rx, ry, rz, vx, vy, vz)
         implicit none
         real(DP), intent(in)  :: mu
            !! Gravitational constant
         real(DP), intent(in)  :: a
            !! semimajor axis
         real(DP), intent(in) :: ie
            !! eccentricity input. If <0 then 0.0 is used
         real(DP), intent(in)  :: inc_deg
            !! inclination (degrees)
         real(DP), intent(in)  :: capom_deg
            !! longitude of ascending node (degrees)
         real(DP), intent(in)  :: omega_deg
            !! argument of periapsis (degrees)
         real(DP), intent(in)  :: capm_deg
            !! mean anomaly (degrees)
         real(DP), intent(out) :: rx, ry, rz
            !! Position vector
         real(DP), intent(out) :: vx, vy, vz
            !! Velocity vector
      end subroutine swiftest_orbel_el2xv

      module subroutine swiftest_orbel_el2xv_vec(self, cb)
         implicit none
         class(swiftest_body), intent(inout) :: self 
            !! Swiftest body object
         class(swiftest_cb),   intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_orbel_el2xv_vec

      pure module subroutine swiftest_orbel_scget(angle, sx, cx)
         !$omp declare simd(swiftest_orbel_scget)
         implicit none
         real(DP), intent(in)  :: angle
         real(DP), intent(out) :: sx, cx
      end subroutine swiftest_orbel_scget

      pure elemental module subroutine swiftest_orbel_xv2aeq(mu, rx, ry, rz, vx, vy, vz, a, e, q)
         !$omp declare simd(swiftest_orbel_xv2aeq)
         implicit none
         real(DP), intent(in)  :: mu       
            !! Gravitational constant
         real(DP), intent(in)  :: rx,ry,rz 
            !! Position vector
         real(DP), intent(in)  :: vx,vy,vz 
            !! Velocity vector
         real(DP), intent(out) :: a        
            !! semimajor axis
         real(DP), intent(out) :: e        
            !! eccentricity
         real(DP), intent(out) :: q        
            !! periapsis
      end subroutine swiftest_orbel_xv2aeq

      pure module subroutine swiftest_orbel_xv2aqt(mu, rx, ry, rz, vx, vy, vz, a, q, capm, tperi)
         !$omp declare simd(swiftest_orbel_xv2aqt)
         implicit none
         real(DP), intent(in)  :: mu       
            !! Gravitational constant
         real(DP), intent(in)  :: rx,ry,rz 
            !! Position vector
         real(DP), intent(in)  :: vx,vy,vz 
            !! Velocity vector
         real(DP), intent(out) :: a        
            !! semimajor axis
         real(DP), intent(out) :: q        
            !! periapsis
         real(DP), intent(out) :: capm     
            !! mean anomaly
         real(DP), intent(out) :: tperi    
            !! time of pericenter passage
      end subroutine swiftest_orbel_xv2aqt

      pure elemental module subroutine swiftest_orbel_xv2el(mu, rx, ry, rz, vx, vy, vz, &
                                                      a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
         implicit none
         real(DP), intent(in)  :: mu    
            !! Gravitational constant
         real(DP), intent(in)  :: rx,ry,rz 
            !! Position vector
         real(DP), intent(in)  :: vx,vy,vz 
            !! Velocity vector
         real(DP), intent(out) :: a     
            !! semimajor axis
         real(DP), intent(out) :: e     
            !! eccentricity
         real(DP), intent(out) :: inc   
            !! inclination
         real(DP), intent(out) :: capom 
            !! longitude of ascending node
         real(DP), intent(out) :: omega 
            !! argument of periapsis
         real(DP), intent(out) :: capm  
            !! mean anomaly
         real(DP), intent(out) :: varpi 
            !! longitude of periapsis
         real(DP), intent(out) :: lam   
            !! mean longitude
         real(DP), intent(out) :: f     
            !! true anomaly
         real(DP), intent(out) :: cape  
            !! eccentric anomaly (eccentric orbits)
         real(DP), intent(out) :: capf  
            !! hyperbolic anomaly (hyperbolic orbits)
      end subroutine swiftest_orbel_xv2el

      module subroutine swiftest_orbel_xv2el_vec(self, cb)
         implicit none
         class(swiftest_body), intent(inout) :: self 
            !! Swiftest body object
         class(swiftest_cb),   intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_orbel_xv2el_vec

      module subroutine swiftest_util_setup_body(self, n, param)
         implicit none
         class(swiftest_body),       intent(inout) :: self  
            !! Swiftest body object
         integer(I4B),               intent(in)    :: n     
            !! Number of particles to allocate space for
         class(swiftest_parameters), intent(in)    :: param 
            !! Current run configuration parameters
      end subroutine swiftest_util_setup_body

      module subroutine swiftest_util_setup_construct_system(nbody_system, param)
         implicit none
         class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system 
            !! Swiftest nbody_system object
         class(swiftest_parameters),                intent(inout) :: param        
            !! Current run configuration parameters
      end subroutine swiftest_util_setup_construct_system

      module subroutine swiftest_util_setup_initialize_particle_info_system(self, param)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self  
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine swiftest_util_setup_initialize_particle_info_system

      module subroutine swiftest_util_setup_initialize_system(self, system_history, param)
         implicit none
         class(swiftest_nbody_system),              intent(inout) :: self           
            !! Swiftest nbody_system object
         class(swiftest_storage),      allocatable, intent(inout) :: system_history 
            !! Stores the system history between output dumps
         class(swiftest_parameters),                intent(inout) :: param          
            !! Current run configuration parameters
      end subroutine swiftest_util_setup_initialize_system

      module subroutine swiftest_util_setup_pl(self, n, param)
         implicit none
         class(swiftest_pl),         intent(inout) :: self  
            !! Swiftest massive body object
         integer(I4B),               intent(in)    :: n     
            !! Number of particles to allocate space for
         class(swiftest_parameters), intent(in)    :: param 
            !! Current run configuration parameters
      end subroutine swiftest_util_setup_pl

      module subroutine swiftest_util_setup_tp(self, n, param)
         implicit none
         class(swiftest_tp),         intent(inout) :: self  
            !! Swiftest test particle object
         integer(I4B),               intent(in)    :: n     
            !! Number of particles to allocate space for
         class(swiftest_parameters), intent(in)    :: param 
            !! Current run configuration parametersr
      end subroutine swiftest_util_setup_tp

      module subroutine swiftest_user_kick_getacch_body(self, nbody_system, param, t, lbeg)
         implicit none
         class(swiftest_body),              intent(inout) :: self   
            !! Swiftest massive body particle data structure
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody_system_object
         class(swiftest_parameters),        intent(inout) :: param  
            !! Current run configuration parameters 
         real(DP),                          intent(in)    :: t      
            !! Current time
         logical,                           intent(in)    :: lbeg   
            !! Optional argument that determines whether or not this is the beginning or end of the step
      end subroutine swiftest_user_kick_getacch_body
   end interface

   interface util_append
      module subroutine swiftest_util_append_arr_info(arr, source, nold, lsource_mask)
         implicit none
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout)        :: arr          
            !! Destination array 
         type(swiftest_particle_info), dimension(:), allocatable, intent(in)           :: source       
            !! Array to append 
         integer(I4B),                                            intent(in), optional :: nold         
            !! Extent of original array. If passed, the source array will begin at arr(nold+1). 
            !! Otherwise, the size of arr will be used.
         logical,                      dimension(:),              intent(in), optional :: lsource_mask 
            !! Logical mask indicating which elements to append to
      end subroutine swiftest_util_append_arr_info

      module subroutine swiftest_util_append_arr_kin(arr, source, nold, lsource_mask)
         implicit none
         type(swiftest_kinship), dimension(:), allocatable, intent(inout)        :: arr          
            !! Destination array 
         type(swiftest_kinship), dimension(:), allocatable, intent(in)           :: source       
            !! Array to append 
         integer(I4B),                                      intent(in), optional :: nold         
            !! Extent of original array. If passed, the source array will begin at arr(nold+1). 
            !! Otherwise, the size of arr will be used.
         logical,                dimension(:),              intent(in), optional :: lsource_mask 
            !! Logical mask indicating which elements to append to
      end subroutine swiftest_util_append_arr_kin
   end interface

   interface
      module subroutine swiftest_util_append_body(self, source, lsource_mask)
         implicit none
         class(swiftest_body),  intent(inout) :: self          
            !! Swiftest body object
         class(swiftest_body),  intent(in)    :: source        
            !! Source object to append
         logical, dimension(:), intent(in)    :: lsource_mask  
            !! Logical mask indicating which elements to append to
      end subroutine swiftest_util_append_body

      module subroutine swiftest_util_append_pl(self, source, lsource_mask)
         implicit none
         class(swiftest_pl),    intent(inout) :: self         
            !! Swiftest massive body object
         class(swiftest_body),  intent(in)    :: source       
            !! Source object to append
         logical, dimension(:), intent(in)    :: lsource_mask 
            !! Logical mask indicating which elements to append to
      end subroutine swiftest_util_append_pl
   
      module subroutine swiftest_util_append_tp(self, source, lsource_mask)
         implicit none
         class(swiftest_tp),    intent(inout) :: self         
            !! Swiftest test particle object
         class(swiftest_body),  intent(in)    :: source       
            !! Source object to append
         logical, dimension(:), intent(in)    :: lsource_mask 
            !! Logical mask indicating which elements to append to
      end subroutine swiftest_util_append_tp

      module subroutine swiftest_util_coord_b2h_pl(self, cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_coord_b2h_pl

      module subroutine swiftest_util_coord_b2h_tp(self, cb)
         implicit none
         class(swiftest_tp), intent(inout) :: self 
            !! Swiftest test particle object
         class(swiftest_cb), intent(in)    :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_coord_b2h_tp

      module subroutine swiftest_util_coord_h2b_pl(self, cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_coord_h2b_pl

      module subroutine swiftest_util_coord_h2b_tp(self, cb)
         implicit none
         class(swiftest_tp), intent(inout) :: self 
            !! Swiftest test particle object
         class(swiftest_cb), intent(in)    :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_coord_h2b_tp

      module subroutine swiftest_util_coord_vb2vh_pl(self, cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_coord_vb2vh_pl
   
      module subroutine swiftest_util_coord_vb2vh_tp(self, vbcb)
         implicit none
         class(swiftest_tp),     intent(inout) :: self 
            !! Swiftest test particle object
         real(DP), dimension(:), intent(in)    :: vbcb 
            !! Barycentric velocity of the central body
      end subroutine swiftest_util_coord_vb2vh_tp
   
      module subroutine swiftest_util_coord_vh2vb_pl(self, cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_coord_vh2vb_pl
   
      module subroutine swiftest_util_coord_vh2vb_tp(self, vbcb)
         implicit none
         class(swiftest_tp),     intent(inout) :: self 
            !! Swiftest test particle object
         real(DP), dimension(:), intent(in)    :: vbcb 
            !! Barycentric velocity of the central body
      end subroutine swiftest_util_coord_vh2vb_tp

      module subroutine swiftest_util_coord_rh2rb_pl(self, cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_coord_rh2rb_pl

      module subroutine swiftest_util_coord_rh2rb_tp(self, cb)
         implicit none
         class(swiftest_tp), intent(inout) :: self 
            !! Swiftest test particle object
         class(swiftest_cb), intent(in) :: cb      
            !! Swiftest central body object
      end subroutine swiftest_util_coord_rh2rb_tp

      module subroutine swiftest_util_copy_particle_info(self, source)
         implicit none
         class(swiftest_particle_info),  intent(inout) :: self
         class(swiftest_particle_info),  intent(in)    :: source
      end subroutine swiftest_util_copy_particle_info

      module subroutine swiftest_util_copy_particle_info_arr(source, dest, idx)
         implicit none
         class(swiftest_particle_info), dimension(:), intent(in)             :: source 
            !! Source object to copy into
         class(swiftest_particle_info), dimension(:), intent(inout)          :: dest   
            !! Swiftest body object with particle metadata information object
         integer(I4B),                  dimension(:), intent(in),   optional :: idx    
            !! Optional array of indices to draw the source object
      end subroutine swiftest_util_copy_particle_info_arr

      module subroutine swiftest_util_dealloc_body(self)
         implicit none
         class(swiftest_body),  intent(inout) :: self
      end subroutine swiftest_util_dealloc_body

      module subroutine swiftest_util_dealloc_kin(self)
         implicit none
         class(swiftest_kinship), intent(inout) :: self 
            !! Swiftest kinship object
      end subroutine swiftest_util_dealloc_kin

      module subroutine swiftest_util_dealloc_cb(self)
         implicit none
         class(swiftest_cb), intent(inout) :: self 
            !! Swiftest central body object
      end subroutine swiftest_util_dealloc_cb

      module subroutine swiftest_util_dealloc_pl(self)
         implicit none
         class(swiftest_pl), intent(inout) :: self
      end subroutine swiftest_util_dealloc_pl

      module subroutine swiftest_util_dealloc_storage(self)
         implicit none
         class(swiftest_storage), intent(inout) :: self 
            !! Swiftest storage object
      end subroutine swiftest_util_dealloc_storage

      module subroutine swiftest_util_dealloc_system(self)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self
      end subroutine swiftest_util_dealloc_system

      module subroutine swiftest_util_dealloc_tp(self)
         implicit none
         class(swiftest_tp), intent(inout) :: self
      end subroutine swiftest_util_dealloc_tp

      module subroutine swiftest_util_fill_body(self, inserts, lfill_list)
         implicit none
         class(swiftest_body),  intent(inout) :: self       
            !! Swiftest body object
         class(swiftest_body),  intent(in)    :: inserts    
            !! Swiftest body object to be inserted
         logical, dimension(:), intent(in)    :: lfill_list 
            !! Logical array of bodies to merge into the keeps
      end subroutine swiftest_util_fill_body

      module subroutine swiftest_util_fill_pl(self, inserts, lfill_list)
         implicit none
         class(swiftest_pl),    intent(inout) :: self       
            !! Swiftest massive body object
         class(swiftest_body),  intent(in)    :: inserts    
            !! Swiftest body object to be inserted
         logical, dimension(:), intent(in)    :: lfill_list 
            !! Logical array of bodies to merge into the keeps
      end subroutine swiftest_util_fill_pl

      module subroutine swiftest_util_fill_tp(self, inserts, lfill_list)
         implicit none
         class(swiftest_tp),    intent(inout) :: self       
            !! Swiftest test particle object
         class(swiftest_body),  intent(in)    :: inserts    
            !! Swiftest body object to be inserted
         logical, dimension(:), intent(in)    :: lfill_list 
            !! Logical array of bodies to merge into the keeps
      end subroutine swiftest_util_fill_tp
   end interface

   interface util_fill
      module subroutine swiftest_util_fill_arr_info(keeps, inserts, lfill_list)
         implicit none
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: keeps      
            !! Array of values to keep 
         type(swiftest_particle_info), dimension(:), allocatable, intent(in)    :: inserts    
            !! Array of values to insert into keep
         logical,                      dimension(:),              intent(in)    :: lfill_list 
            !! Logical array of bodies to merge into the keeps
      end subroutine swiftest_util_fill_arr_info

      module subroutine swiftest_util_fill_arr_kin(keeps, inserts, lfill_list)
         implicit none
         type(swiftest_kinship), dimension(:), allocatable, intent(inout) :: keeps      
            !! Array of values to keep 
         type(swiftest_kinship), dimension(:), allocatable, intent(in)    :: inserts    
            !! Array of values to insert into keep
         logical,                dimension(:),              intent(in)    :: lfill_list 
            !! Logical array of bodies to merge into the keeps
      end subroutine swiftest_util_fill_arr_kin
   end interface

   interface
      pure module subroutine swiftest_util_flatten_eucl_ij_to_k(n, i, j, k)
         !$omp declare simd(swiftest_util_flatten_eucl_ij_to_k)
         implicit none
         integer(I4B), intent(in)  :: n 
            !! Number of bodies
         integer(I4B), intent(in)  :: i 
            !! Index of the ith body
         integer(I4B), intent(in)  :: j 
            !! Index of the jth body
         integer(I8B), intent(out) :: k 
            !! Index of the flattened matrix
      end subroutine swiftest_util_flatten_eucl_ij_to_k

      pure module subroutine swiftest_util_flatten_eucl_k_to_ij(n, k, i, j)
         implicit none
         integer(I4B), intent(in)  :: n 
            !! Number of bodies
         integer(I8B), intent(in)  :: k 
            !! Index of the flattened matrix
         integer(I4B), intent(out) :: i 
            !! Index of the ith body
         integer(I4B), intent(out) :: j 
            !! Index of the jth body
      end subroutine swiftest_util_flatten_eucl_k_to_ij

      module subroutine swiftest_util_flatten_eucl_plpl(self, param)
         implicit none
         class(swiftest_pl),         intent(inout) :: self  
            !! Swiftest massive body object
         class(swiftest_parameters), intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine

      module subroutine swiftest_util_flatten_eucl_pltp(self, pl, param)
         implicit none
         class(swiftest_tp),         intent(inout) :: self  
            !! Swiftest test particle object
         class(swiftest_pl),         intent(in)    :: pl    
            !! Swiftest massive body object
         class(swiftest_parameters), intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine

      module subroutine swiftest_util_get_energy_and_momentum_system(self, param)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self     
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(in)    :: param    
            !! Current run configuration parameters
      end subroutine swiftest_util_get_energy_and_momentum_system

      module subroutine swiftest_util_get_idvalues_system(self, idvals)
         implicit none
         class(swiftest_nbody_system),       intent(in)  :: self   
            !! Encounter snapshot object
         integer(I4B), dimension(:), allocatable, intent(out) :: idvals 
            !! Array of all id values saved in this snapshot
      end subroutine swiftest_util_get_idvalues_system
   end interface

   interface swiftest_util_get_potential_energy
      module subroutine swiftest_util_get_potential_energy_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, rb, pe)
         implicit none
         integer(I4B),                 intent(in)  :: npl
         integer(I8B),                 intent(in)  :: nplpl
         integer(I4B), dimension(:,:), intent(in)  :: k_plpl
         logical,      dimension(:),   intent(in)  :: lmask
         real(DP),                     intent(in)  :: GMcb
         real(DP),     dimension(:),   intent(in)  :: Gmass
         real(DP),     dimension(:),   intent(in)  :: mass
         real(DP),     dimension(:,:), intent(in)  :: rb
         real(DP),                     intent(out) :: pe
      end subroutine swiftest_util_get_potential_energy_flat
   
      module subroutine swiftest_util_get_potential_energy_triangular(npl, lmask, GMcb, Gmass, mass, rb, pe)
         implicit none
         integer(I4B),                 intent(in)  :: npl
         logical,      dimension(:),   intent(in)  :: lmask
         real(DP),                     intent(in)  :: GMcb
         real(DP),     dimension(:),   intent(in)  :: Gmass
         real(DP),     dimension(:),   intent(in)  :: mass
         real(DP),     dimension(:,:), intent(in)  :: rb
         real(DP),                     intent(out) :: pe
      end subroutine swiftest_util_get_potential_energy_triangular
   end interface

   interface
      module subroutine swiftest_util_get_vals_storage(self, idvals, tvals)
         class(swiftest_storage),              intent(in)  :: self   
            !! Swiftest storage object
         integer(I4B), dimension(:), allocatable, intent(out) :: idvals 
            !! Array of all id values in all snapshots
         real(DP),     dimension(:), allocatable, intent(out) :: tvals  
            !! Array of all time values in all snapshots
      end subroutine swiftest_util_get_vals_storage

      module subroutine swiftest_util_index_array(ind_arr, n)
         implicit none
         integer(I4B), dimension(:), allocatable, intent(inout) :: ind_arr 
            !! Index array. Input is a pre-existing index array where n /= size(ind_arr). 
            !! Output is a new index array ind_arr = [1, 2, ... n]
         integer(I4B),                            intent(in)    :: n       
            !! The new size of the index array
      end subroutine swiftest_util_index_array

      module subroutine swiftest_util_index_map_storage(self)
         implicit none
         class(swiftest_storage), intent(inout) :: self 
            !! Swiftest storage object
      end subroutine swiftest_util_index_map_storage

      module subroutine swiftest_util_make_impactors_pl(self, idx)
         implicit none
         class(swiftest_pl),         intent(inout) :: self  
            !! Massive body object
         integer(I4B), dimension(:), intent(in)    :: idx 
            !! Array holding the indices of the two bodies involved in the collision)
      end subroutine swiftest_util_make_impactors_pl

      module subroutine swiftest_util_peri(n,m, r, v, atp, q, isperi)
         implicit none
         integer(I4B),                 intent(in)    :: n      
            !! Number of bodies
         real(DP),     dimension(:),   intent(in)    :: m      
            !! Mass term (mu for HELIO coordinates, and Gmtot for BARY)
         real(DP),     dimension(:,:), intent(in)    :: r      
            !! Position vectors (rh for HELIO coordinates, rb for BARY)
         real(DP),     dimension(:,:), intent(in)    :: v      
            !! Position vectors (vh for HELIO coordinates, rb for BARY)
         real(DP),     dimension(:),   intent(out)   :: atp    
            !! Semimajor axis 
         real(DP),     dimension(:),   intent(out)   :: q      
            !! Periapsis
         integer(I4B), dimension(:),   intent(inout) :: isperi 
            !! Periapsis passage flag
      end subroutine swiftest_util_peri

      module subroutine swiftest_util_peri_body(self, nbody_system, param)
         implicit none
         class(swiftest_body),         intent(inout) :: self   
            !! SyMBA massive body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),   intent(in)    :: param  
            !! Current run configuration parameters
      end subroutine swiftest_util_peri_body

      module subroutine swiftest_util_peri_tp(self, nbody_system, param) 
         implicit none
         class(swiftest_tp),                intent(inout) :: self   
            !! Swiftest test particle object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(in)    :: param  
            !! Current run configuration parameters
      end subroutine swiftest_util_peri_tp

      module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
         implicit none
         class(swiftest_pl),           intent(inout) :: self   
            !! Swiftest massive body object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),   intent(inout) :: param  
            !! Current run configuration parameters 
      end subroutine swiftest_util_rearray_pl

      module subroutine swiftest_util_rearray_tp(self, nbody_system, param)
         implicit none
         class(swiftest_tp),           intent(inout) :: self   
            !! Swiftest test particle object
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(swiftest_parameters),   intent(inout) :: param  
            !! Current run configuration parameters 
      end subroutine swiftest_util_rearray_tp

      module subroutine swiftest_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 swiftest_util_rescale_system

      module subroutine swiftest_util_reset_kinship_pl(self, idx)
         implicit none
         class(swiftest_pl),         intent(inout) :: self 
            !! SyMBA massive body object
         integer(I4B), dimension(:), intent(in)    :: idx  
            !! Index array of bodies to reset
      end subroutine swiftest_util_reset_kinship_pl

   end interface


   interface util_resize
      module subroutine swiftest_util_resize_arr_info(arr, nnew)
         implicit none
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr  
            !! Array to resize
         integer(I4B),                                            intent(in)    :: nnew 
            !! New size
      end subroutine swiftest_util_resize_arr_info

      module subroutine swiftest_util_resize_arr_kin(arr, nnew)
         implicit none
         type(swiftest_kinship), dimension(:), allocatable, intent(inout) :: arr  
            !! Array to resize
         integer(I4B),                                       intent(in)    :: nnew 
            !! New size
      end subroutine swiftest_util_resize_arr_kin
   end interface

   interface
      module subroutine swiftest_util_resize_body(self, nnew)
         implicit none
         class(swiftest_body), intent(inout) :: self 
            !! Swiftest body object
         integer(I4B),         intent(in)    :: nnew 
            !! New size neded
      end subroutine swiftest_util_resize_body

      module subroutine swiftest_util_resize_pl(self, nnew)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         integer(I4B),       intent(in)    :: nnew 
            !! New size neded
      end subroutine swiftest_util_resize_pl

      module subroutine swiftest_util_resize_storage(storage, nold, nnew)
         use base, only : base_storage
         implicit none
         class(base_storage), allocatable, intent(inout) :: storage 
            !! Original storage object
         integer(I4B),                        intent(in)    :: nold    
            !! Old size
         integer(I4B),                        intent(in)    :: nnew    
            !! New size
      end subroutine swiftest_util_resize_storage 

      module subroutine swiftest_util_resize_tp(self, nnew)
         implicit none
         class(swiftest_tp), intent(inout) :: self 
            !! Swiftest test particle object
         integer(I4B),       intent(in)    :: nnew 
            !! New size neded
      end subroutine swiftest_util_resize_tp

      module subroutine swiftest_util_save_discard_body(self, ldiscard_mask, nbody_system, snapshot)
         implicit none
         class(swiftest_body), intent(inout) :: self 
            !! Swiftest body object
         logical, dimension(:), intent(in) :: ldiscard_mask 
            !! Logical mask indicating which elements to discard
         class(swiftest_nbody_system), intent(inout) :: nbody_system 
            !! Swiftest nbody system object
         class(base_nbody_system), intent(inout) :: snapshot
            !! Nbody system object to store the snapshot of the discard system
      end subroutine swiftest_util_save_discard_body

      module subroutine swiftest_util_set_beg_end_pl(self, rbeg, rend, vbeg)
         implicit none
         class(swiftest_pl),       intent(inout)          :: self 
            !! Swiftest massive body object
         real(DP), dimension(:,:), intent(in),   optional :: rbeg 
            !! Position vectors at beginning of step
         real(DP), dimension(:,:), intent(in),   optional :: rend 
            !! Positions vectors at end of step
         real(DP), dimension(:,:), intent(in),   optional :: vbeg 
            !! vbeg is an unused variable to keep this method forward compatible with RMVS
      end subroutine swiftest_util_set_beg_end_pl

      module subroutine swiftest_util_set_ir3h(self)
         implicit none
         class(swiftest_body), intent(inout) :: self 
            !! Swiftest body object
      end subroutine swiftest_util_set_ir3h

      module subroutine swiftest_util_set_msys(self)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self 
            !! Swiftest nbody_system object
      end subroutine swiftest_util_set_msys

      module subroutine swiftest_util_set_mu_pl(self, cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_set_mu_pl

      module subroutine swiftest_util_set_mu_tp(self, cb)
         implicit none
         class(swiftest_tp), intent(inout) :: self 
            !! Swiftest test particle object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_set_mu_tp

      module subroutine swiftest_util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, &
                                               origin_rh, origin_vh, discard_time, discard_rh, discard_vh, discard_body_id)
         implicit none
         class(swiftest_particle_info), intent(inout)           :: self
         character(len=*),              intent(in),    optional :: name            
            !! Non-unique name
         character(len=*),              intent(in),    optional :: particle_type   
            !! String containing a description of the particle type (e.g. Central Body, Massive Body, Test Particle)
         character(len=*),              intent(in),    optional :: status          
            !! Particle status description: Active, Merged, Fragmented, etc.
         character(len=*),              intent(in),    optional :: origin_type     
            !! String containing a description of the origin of the particle (e.g. Initial Conditions, Disruption, etc.)
         real(DP),                      intent(in),    optional :: origin_time     
            !! The time of the particle's formation
         integer(I4B),                  intent(in),    optional :: collision_id    
            !! The ID fo the collision that formed the particle
         real(DP), dimension(:),        intent(in),    optional :: origin_rh       
            !! The heliocentric distance vector at the time of the particle's formation
         real(DP), dimension(:),        intent(in),    optional :: origin_vh       
            !! The heliocentric velocity vector at the time of the particle's formation
         real(DP),                      intent(in),    optional :: discard_time    
            !! The time of the particle's discard
         real(DP), dimension(:),        intent(in),    optional :: discard_rh      
            !! The heliocentric distance vector at the time of the particle's discard
         real(DP), dimension(:),        intent(in),    optional :: discard_vh      
            !! The heliocentric velocity vector at the time of the particle's discard
         integer(I4B),                  intent(in),    optional :: discard_body_id 
            !! The id of the other body involved in the discard (0 if no other body involved)
      end subroutine swiftest_util_set_particle_info

      module subroutine swiftest_util_set_rhill(self,cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_set_rhill

      module subroutine swiftest_util_set_renc_I4B(self, scale)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         integer(I4B),       intent(in)    :: scale 
            !! Input scale factor (multiplier of Hill's sphere size)
      end subroutine swiftest_util_set_renc_I4B

      module subroutine swiftest_util_set_renc_DP(self, scale)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         real(DP),           intent(in)    :: scale 
            !! Input scale factor (multiplier of Hill's sphere size)
      end subroutine swiftest_util_set_renc_DP

      module subroutine swiftest_util_set_rhill_approximate(self,cb)
         implicit none
         class(swiftest_pl), intent(inout) :: self 
            !! Swiftest massive body object
         class(swiftest_cb), intent(inout) :: cb   
            !! Swiftest central body object
      end subroutine swiftest_util_set_rhill_approximate

      module subroutine swiftest_util_snapshot_save(storage, snapshot)
         use base, only : base_storage
         implicit none
         class(base_storage), allocatable, intent(inout) :: storage  
            !! Storage ncounter storage object
         class(*),                         intent(in)    :: snapshot 
            !! Object to snapshot
      end subroutine swiftest_util_snapshot_save

      module subroutine swiftest_util_snapshot_system(self, param, nbody_system, t, arg)
         implicit none
         class(swiftest_storage),      intent(inout)        :: self            
            !! Swiftest storage object
         class(swiftest_parameters),   intent(inout)        :: param           
            !! Current run configuration parameters
         class(swiftest_nbody_system), intent(inout)        :: nbody_system    
            !! Swiftest nbody system object to store
         real(DP),                     intent(in), optional :: t               
            !! Time of snapshot if different from nbody_system time
         character(*),                 intent(in), optional :: arg             
            !! Optional argument (needed for extended storage type used in collision snapshots)
      end subroutine swiftest_util_snapshot_system
   end interface


   interface util_sort_rearrange
      module subroutine swiftest_util_sort_rearrange_arr_info(arr, ind, n)
         implicit none
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr 
            !! Destination array 
         integer(I4B),                 dimension(:),              intent(in)    :: ind 
            !! Index to rearrange against
         integer(I4B),                                            intent(in)    :: n   
            !! Number of elements in arr and ind to rearrange
      end subroutine swiftest_util_sort_rearrange_arr_info

      pure module subroutine swiftest_util_sort_rearrange_arr_kin(arr, ind, n)
         implicit none
         type(swiftest_kinship), dimension(:), allocatable, intent(inout) :: arr 
            !! Destination array 
         integer(I4B),           dimension(:),              intent(in)    :: ind 
            !! Index to rearrange against
         integer(I4B),                                      intent(in)    :: n   
            !! Number of elements in arr and ind to rearrange
      end subroutine swiftest_util_sort_rearrange_arr_kin

   end interface util_sort_rearrange

   interface
      module subroutine swiftest_util_sort_rearrange_body(self, ind)
         implicit none
         class(swiftest_body),               intent(inout) :: self 
            !! Swiftest body object
         integer(I4B),         dimension(:), intent(in)    :: ind  
            !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
      end subroutine swiftest_util_sort_rearrange_body

      module subroutine swiftest_util_sort_rearrange_pl(self, ind)
         implicit none
         class(swiftest_pl),                 intent(inout) :: self 
            !! Swiftest massive body object
         integer(I4B),         dimension(:), intent(in)    :: ind  
            !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
      end subroutine swiftest_util_sort_rearrange_pl

      module subroutine swiftest_util_sort_rearrange_tp(self, ind)
         implicit none
         class(swiftest_tp),                 intent(inout) :: self 
            !! Swiftest test particle object
         integer(I4B),         dimension(:), intent(in)    :: ind  
            !! Index array used to restructure the body (should contain all 1:n index values in the desired order)
      end subroutine swiftest_util_sort_rearrange_tp

      module subroutine swiftest_util_sort_body(self, sortby, ascending)
         implicit none
         class(swiftest_body), intent(inout) :: self      
            !! Swiftest body object
         character(*),         intent(in)    :: sortby    
            !! Sorting attribute
         logical,              intent(in)    :: ascending 
            !! Logical flag indicating whether or not the sorting should be in ascending or descending order
      end subroutine swiftest_util_sort_body

      module subroutine swiftest_util_sort_pl(self, sortby, ascending)
         implicit none
         class(swiftest_pl), intent(inout) :: self      
            !! Swiftest body object
         character(*),       intent(in)    :: sortby    
            !! Sorting attribute
         logical,            intent(in)    :: ascending 
            !! Logical flag indicating whether or not the sorting should be in ascending or descending order
      end subroutine swiftest_util_sort_pl

      module subroutine swiftest_util_sort_tp(self, sortby, ascending)
         implicit none
         class(swiftest_tp), intent(inout) :: self      
            !! Swiftest body object
         character(*),       intent(in)    :: sortby    
            !! Sorting attribute
         logical,            intent(in)    :: ascending 
            !! Logical flag indicating whether or not the sorting should be in ascending or descending order
      end subroutine swiftest_util_sort_tp

   end interface

   interface util_spill
      module subroutine swiftest_util_spill_arr_info(keeps, discards, lspill_list, ldestructive)
         implicit none
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: keeps        
            !! Array of values to keep 
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: discards     
            !! Array of discards
         logical,                      dimension(:),              intent(in)    :: lspill_list  
            !! Logical array of bodies to spill into the discardss
         logical,                                                 intent(in)    :: ldestructive 
            !! Logical flag indicating whether or not this operation should alter the keeps array or not
      end subroutine swiftest_util_spill_arr_info

      module subroutine swiftest_util_spill_arr_kin(keeps, discards, lspill_list, ldestructive)
         implicit none
         type(swiftest_kinship), dimension(:), allocatable, intent(inout) :: keeps        
            !! Array of values to keep 
         type(swiftest_kinship), dimension(:), allocatable, intent(inout) :: discards     
            !! Array of discards
         logical,                dimension(:),              intent(in)    :: lspill_list  
            !! Logical array of bodies to spill into the discardss
         logical,                                           intent(in)    :: ldestructive 
            !! Logical flag indicating whether or not this operation should alter the keeps array or not
      end subroutine swiftest_util_spill_arr_kin
   end interface

   interface 
      module subroutine swiftest_util_spill_body(self, discards, lspill_list, ldestructive)
         implicit none
         class(swiftest_body),  intent(inout) :: self         
            !! Swiftest body object
         class(swiftest_body),  intent(inout) :: discards     
            !! Discarded object 
         logical, dimension(:), intent(in)    :: lspill_list  
            !! Logical array of bodies to spill into the discards
         logical,               intent(in)    :: ldestructive 
            !! Logical flag indicating whether or not this operation should alter the keeps array or not
      end subroutine swiftest_util_spill_body

      module subroutine swiftest_util_spill_pl(self, discards, lspill_list, ldestructive)
         implicit none
         class(swiftest_pl),    intent(inout) :: self         
            !! Swiftest massive body object
         class(swiftest_body),  intent(inout) :: discards     
            !! Discarded object 
         logical, dimension(:), intent(in)    :: lspill_list  
            !! Logical array of bodies to spill into the discards
         logical,               intent(in)    :: ldestructive 
            !! Logical flag indicating whether or not this operation should alter the keeps array or not
      end subroutine swiftest_util_spill_pl

      module subroutine swiftest_util_spill_tp(self, discards, lspill_list, ldestructive)
         implicit none
         class(swiftest_tp),    intent(inout) :: self         
            !! Swiftest test particle object
         class(swiftest_body),  intent(inout) :: discards     
            !! Discarded object 
         logical, dimension(:), intent(in)    :: lspill_list  
            !! Logical array of bodies to spill into the discards
         logical,               intent(in)    :: ldestructive 
            !! Logical flag indicating whether or not this operation should alter the keeps array or not
      end subroutine swiftest_util_spill_tp

   end interface

   interface
      module subroutine swiftest_util_valid_id_system(self, param)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self  
            !! Swiftest nbody system object
         class(swiftest_parameters),        intent(inout) :: param 
            !! Current run configuration parameters
      end subroutine swiftest_util_valid_id_system

      module subroutine swiftest_util_version()
         implicit none
      end subroutine swiftest_util_version
   end interface

#ifdef COARRAY
   interface
      module subroutine swiftest_coarray_balance_system(self, param)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self
            !! Swiftest nbody system 
         class(swiftest_parameters),   intent(inout) :: param        
            !! Current run configuration parameters 
      end subroutine swiftest_coarray_balance_system

      module subroutine swiftest_coarray_collect_system(self, param)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self
            !! Swiftest nbody system 
         class(swiftest_parameters),   intent(inout) :: param        
            !! Current run configuration parameters 
      end subroutine swiftest_coarray_collect_system

      module subroutine swiftest_coarray_distribute_system(self, param)
         implicit none
         class(swiftest_nbody_system), intent(inout) :: self
            !! Swiftest nbody system 
         class(swiftest_parameters),   intent(inout) :: param        
            !! Current run configuration parameters 
      end subroutine swiftest_coarray_distribute_system
   end interface

   interface coclone
      module subroutine swiftest_coarray_component_clone_info(var,src_img)
         implicit none
         type(swiftest_particle_info), intent(inout) :: var
         integer(I4B), intent(in),optional :: src_img
      end subroutine swiftest_coarray_component_clone_info

      module subroutine swiftest_coarray_component_clone_info_arr1D(var,src_img)
         implicit none
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: var
         integer(I4B), intent(in),optional :: src_img
      end subroutine swiftest_coarray_component_clone_info_arr1D

      module subroutine swiftest_coarray_coclone_param(param)
         implicit none
         type(swiftest_parameters),intent(inout),codimension[*]  :: param  
      end subroutine swiftest_coarray_coclone_param
   end interface

   interface cocollect
      module subroutine swiftest_coarray_component_collect_info_arr1D(var,dest_img)
         implicit none
         type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: var
         integer(I4B), intent(in),optional :: dest_img
      end subroutine
   end interface

#endif

   contains
      subroutine swiftest_final_kin(self)
         !! author: David A. Minton
         !!
         !! Finalize the swiftest kinship object - deallocates all allocatables
         implicit none
         ! Argument
         type(swiftest_kinship),  intent(inout) :: self 
            !! SyMBA kinship object
         call self%dealloc()
         return
      end subroutine swiftest_final_kin


      subroutine swiftest_final_storage(self)
         !! author: David A. Minton
         !!
         !! Finalizer for the storage data type
         implicit none
         ! Arguments
         type(swiftest_storage) :: self
            !! Swiftest storage object
         call self%dealloc()
         return
      end subroutine swiftest_final_storage

end module swiftest
