C API Reference

This section provides detailed documentation for the C components of NSphere, generated from source code comments using Doxygen and Breathe.

NSphere Core (nsphere.c)

Defines

MYBINIO_H

CLEAN_EXIT(code)

Thread-safe exit macro that properly cleans up allocated resources.

This macro ensures proper resource cleanup before program termination:

Uses OpenMP critical section to ensure only one thread performs cleanup.
Calls cleanup_all_particle_data() and free_local_snap_arrays() to free dynamically allocated memory.
Exits with the specified error code.

Example

if (error_condition) {
    CLEAN_EXIT(1);
}

Parameters:

code – The exit code for the program.

imin(a, b): GLOBAL CONFIGURATION AND MACROS

IF_DEBUG

Conditional compilation macros for feature flags.

These macros provide a cleaner syntax for conditional code blocks that depend on the global feature flags.

Macro for code blocks executed only if g_doDebug is true.

IF_DYNPSI: Macro for code blocks executed only if g_doDynPsi is true.

IF_DYNRANK: Macro for code blocks executed only if g_doDynRank is true.

IF_ALL_PART: Macro for code blocks executed only if g_doAllParticleData is true.

PI: Mathematical constant Pi.

PHYSICAL CONSTANTS AND ASTROPHYSICAL PARAMETERS

Core constants and unit conversion factors for astrophysical calculations

G_CONST: Newton’s gravitational constant in kpc (km/sec)^2/Msun.

RC: Core radius in kpc.

sqr(x)

cube(x)

kmsec_to_kpcmyr: Conversion factor: km/s to kpc/Myr.

HALO_MASS: Default halo mass in solar masses (Msun).

VEL_CONV_SQ: Velocity conversion squared (kpc/Myr)^2 per (km/s)^2.

CLEAN_LOCAL_EXIT(code)

Non-thread-safe exit macro for cleanup and termination.

Performs cleanup and exits in single-threaded contexts:

Frees local snapshot arrays via free_local_snap_arrays().
Cleans up all global particle data arrays via cleanup_all_particle_data().
Exits with the specified error code.

Example

if (error_condition) {
    CLEAN_LOCAL_EXIT(1);
}

Warning

Not thread-safe. Use only in single-threaded contexts or where thread safety is guaranteed by other means.

Parameters:

code – The exit code for the program.

DEBUG_PARTICLE_ID: Energy Debugging and Validation Subsystem

DEBUG_MAX_STEPS

Functions

static int __attribute__ ((unused))

void get_suffixed_filename(const char *base_filename, int with_suffix, char *buffer, size_t bufsize)

Applies the global suffix to a filename.

For .dat files, inserts suffix before the extension. For other files, appends suffix to the end of filename.

void log_message(const char *level, const char *format, ...)

void format_file_size(long size_in_bytes, char *buffer, size_t buffer_size): Formats a byte count into a human-readable string with appropriate units.

int fprintf_bin(FILE *fp, const char *format, ...)

Binary file I/O function declarations.

Writes binary data to a file using a printf-like format string.

These functions provide platform-independent binary file I/O operations similar to fprintf and fscanf, but ensure consistent binary format regardless of host architecture.

BINARY FILE I/O UTILITIES

Parses a format string containing simplified specifiers (d, f, g, e). Writes corresponding arguments from the variadic list (...) as binary data. Integer types (d) are written as int. Floating-point types (f, g, e) are read as double from args but written as float to the file for storage efficiency.

Parameters

fp : FILE* File pointer to write to (must be opened in binary mode). format : const char* Format string with specifiers (d, f, g, e). Other characters are ignored. … : Variable arguments matching the format specifiers.

Returns

int The number of items successfully written according to the format string.

See also

fscanf_bin

Note

Skips optional width/precision specifiers in the format string.

int fscanf_bin(FILE *fp, const char *format, ...)

Reads binary data from a file using a scanf-like format string.

Parses a format string containing simplified specifiers (d, f, g, e). Reads corresponding binary data from the file and stores it in the pointer arguments provided in the variadic list (...). Integer types (d) are read as int. Floating-point types (f, g, e) are read as float from the file but stored into double* arguments provided by the caller.

void cleanup_all_particle_data(void)

Frees all global arrays used for particle data processing.

Frees L_block, Rank_block, R_block, Vrad_block, and chosen. Does NOT free lowestL_* arrays as they are handled elsewhere.

static void free_local_snap_arrays(void)

Frees local arrays used for snapshot processing.

Frees Rank_partdata_snap, R_partdata_snap, Vrad_partdata_snap, L_partdata_snap.

static inline double gravitational_force(double r, int current_rank, int npts, double G_value, double halo_mass_value): Calculates gravitational acceleration at a given radius.

static inline double effective_angular_force(double r, double ell): Computes the effective centrifugal acceleration due to angular momentum.

static inline double gravitational_force_rho_v(double rho, int current_rank, int npts, double G_value, double halo_mass_value)

Alternative gravitational force calculation using transformed coordinates.

Used in the Levi-Civita regularization scheme. Calculates F/m in rho coordinates.

static inline double effective_angular_force_rho_v(double rho, double ell)

Computes effective centrifugal acceleration in transformed coordinates.

Used in the Levi-Civita regularization scheme. Calculates L^2/r^3 in rho coordinates.

static void store_debug_approxE(int snapIndex, double E_value, double time_val): Records the theoretical model energy for a debug snapshot.

static void store_debug_dynE_components(int snapIndex, double totalE, double kinE, double potE, double time_val, double radius_val): Records the actual dynamical energy and its components for a debug snapshot.

static void finalize_debug_energy_output(void)

Writes the collected debug energy comparison data to a file.

This function is called at the end of the simulation if debug mode (g_doDebug) is enabled. It iterates through the stored debug snapshots collected during the simulation and post-processing, writing the following data for the tracked particle (DEBUG_PARTICLE_ID) to data/debug_energy_compare.dat (with the appropriate file suffix applied):

Snapshot index
Simulation time (Myr)
Particle radius (kpc)
Approximate (theoretical) energy (per unit mass)
Dynamic (calculated) energy (per unit mass)
Difference between dynamic and approximate energy
Kinetic energy component (per unit mass)
Potential energy component (per unit mass)

See also

store_debug_approxE

See also

store_debug_dynE_components

Note

The file writing operation is skipped if the simulation is in restart mode and file writes are disabled (skip_file_writes is true).

int double_cmp(const void *a, const void *b)

Comparison function for sorting double values (ascending order).

Used with qsort and other standard library sorting functions.

int compare_particles(const void *a, const void *b)

Comparison function for qsort to order particles based on radius.

Compares two particle entries for sorting based on the first column value.

PARTICLE SORTING AND ORDERING OPERATIONS

Parameters

a : const void* Pointer to the first particle data array (double**). b : const void* Pointer to the second particle data array (double**).

Returns

int Integer less than, equal to, or greater than zero if the radius (first element) of particle a is found, respectively, to be less than, to match, or be greater than the radius of particle b.

Used as a comparison function for qsort, quadsort, and other sorting algorithms. Particles with smaller values in column 0 will be sorted before those with larger values.

Note

Compares based on the first element (index 0) of the column data, which represents the radius.

Note

Assumes the particle data is structured as [particle_index][component_index] when passed via columns array in sorting functions, and compares columns[i][0].

Parameters:

a – Pointer to the first particle entry (as void pointer, expected double**).
b – Pointer to the second particle entry (as void pointer, expected double**).

Returns:

-1 if a<b, 1 if a>b, 0 if equal based on the first column value (radius).

void sort_particles(double **particles, int npts)

Sorts particle data using the default sorting algorithm (g_defaultSortAlg).

Convenience wrapper function for sorting particles with the default algorithm.

void sort_particles_with_alg(double **particles, int npts, const char *sortAlg)

Sorts particle data using the specified sorting algorithm.

Main function for sorting particle data using a specified algorithm.

void sort_rr_psi_arrays(double *rrA_spline, double *psiAarr_spline, int npts)

Sorts radius and psi arrays in tandem based on radius values.

Sorts radius and potential arrays simultaneously to maintain correspondence.

Ensures that the correspondence between radius and potential values is maintained after sorting, which is necessary for spline creation.

double massintegrand(double r, void *params): Integrand for calculating enclosed mass: r^2 * rho(r).

PHYSICS CALCULATION FUNCTIONS

Used in GSL integration routines to compute M(r).

Parameters

r : double Radial coordinate (kpc). params : void* Integration parameters (unused here, but required by GSL).

Returns

double Value of the mass integrand r^2 * rho(r) at radius r.

Note

Assumes a density profile rho(r) proportional to 1 / (1 + (r/RC)^2)^3.

double drhodr(double r)

Calculates the derivative of the density profile with respect to radius.

Calculates the derivative of density with respect to radius.

double Psiintegrand(double rp, void *params)

Integrand for calculating gravitational potential: massintegrand(r) / r.

Calculates the gravitational potential integrand for a given radius.

Used in GSL integration routines to compute Psi(r) via integral from r to infinity.

double evaluatespline(gsl_spline *spline, gsl_interp_accel *acc, double value)

Evaluates a GSL spline at a given value with bounds checking.

Safely evaluates a GSL spline at a given value with robust bounds checking.

Safely evaluates the spline, clamping the input value to the spline’s defined range if necessary to prevent GSL errors.

double fEintegrand(double t, void *params)

Integrand for energy distribution function f(E) calculation using Eddington’s formula.

Calculates the integrand for the distribution function calculation.

Computes 2 * drho/dpsi as a function of t = sqrt(E - Psi_min). Used in GSL integration to find the inner integral in Eddington’s formula.

int cmp_LAI(const void *a, const void *b)

Comparison function for sorting LAndIndex structures by the L member (ascending).

Comparison function for sorting LAndIndex structures by angular momentum.

static int isInteger(const char *str)

Checks if a string represents a valid integer.

Allows an optional leading ‘+’ or ‘-’ sign. Checks if all subsequent characters are digits. Returns false for empty strings or strings containing only a sign.

static int isFloat(const char *str)

Checks if a string is a valid floating point number.

Parameters:: str – The string to check
Returns:: 1 if string is a valid float, 0 otherwise

static void printUsage(const char *prog): Displays detailed usage information for command-line arguments.

COMMAND LINE ARGUMENT PROCESSING

Prints a comprehensive help message to stderr showing all available command-line options, their default values, and brief descriptions. Includes information about integration methods, sorting algorithms, data saving modes, and basic usage examples.

Parameters

prog : const char* The program name to display in the usage message (typically argv[0]).

Returns

None (prints to stderr).

Note

Called when the user specifies —help or when errors occur during argument parsing.

static void errorAndExit(const char *msg, const char *arg, const char *prog)

Displays an error message, suggests —help, and terminates the program.

Formats and prints an error message to stderr, includes a suggestion to use the —help flag for usage information, performs necessary cleanup of allocated resources, and then exits with a non-zero status.

void sort_by_rad(struct PartData *array, int npts)

Sorts particle data in ascending order by radial position.

Sorts an array of PartData structures by their radius value.

Performs a quick sort of the particle data structure array, using the radial position as the primary sorting key. Preserves the original indices to allow tracking particles across time steps.

Includes safety checks for NULL array and invalid counts before proceeding with the sort operation.

Parameters:

array – Array of particle data structures to be sorted
npts – Number of particles in the array
array – Array of PartData structures to sort
npts – Number of elements in the array

static void append_all_particle_data_chunk_to_file(const char *filename, int npts, int block_size, float *L_block, int *Rank_block, float *R_block, float *Vrad_block)

Appends a block of particle data to an output file.

Writes particle data (rank, radius, velocity, angular momentum) in step-major order, where all particles for a single time step are stored contiguously. This layout enables efficient seeking to specific snapshots during post-processing and visualization phases.

Parameters:

filename – Path to the output file
npts – Number of particles
block_size – Number of time steps in this block
L_block – Angular momentum data block (npts * block_size)
Rank_block – Particle rank data block (npts * block_size)
R_block – Radial position data block (npts * block_size)
Vrad_block – Radial velocity data block (npts * block_size)

static void retrieve_all_particle_snapshot(const char *filename, int snap, int npts, int block_size, float *L_out, int *Rank_out, float *R_out, float *Vrad_out)

Retrieves particle data for a specific snapshot from a binary file.

This function reads particle data (rank, radius, velocity, angular momentum) for a specific snapshot number from a binary file (all_particle_data.dat). It allocates temporary thread-local buffers to avoid race conditions when called in parallel. File I/O (seeking and reading) is performed within an OpenMP critical section (file_access) to ensure thread safety. After reading, the data is copied from the local buffers to the output arrays provided by the caller.

Parameters:

filename – Path to the binary data file
snap – Snapshot number to retrieve
npts – Number of particles
block_size – Number of snapshots per block in the file
L_out – Output buffer for angular momentum values
Rank_out – Output buffer for particle ranks
R_out – Output buffer for radial positions
Vrad_out – Output buffer for radial velocities

static int adjust_ntimesteps(int Ntimes_initial, int nout, int dtwrite)

Adjusts the number of timesteps to ensure proper alignment with output intervals.

This function calculates an adjusted number of timesteps (N) such that it is at least the requested initial value (N) and satisfies the constraint: (N - 1) = (M - 1) * k * p for some integer k >= 1, where M (nout) is the number of snapshots, and p (dtwrite) is the write interval.

This ensures that the simulation will produce exactly M output snapshots at intervals of approximately k*p timesteps.

Parameters:

Ntimes_initial – Initially requested number of timesteps (N)
nout – Number of output snapshots needed (M)
dtwrite – Write interval in timesteps (p)

Returns:

Adjusted number of timesteps (N’) satisfying the constraint

static void doMicroLeapfrog(int i, int npts, double r_in, double v_in, double ell, double h, int N, int subSteps, double grav, double *r_out, double *v_out)

Performs a micro-subdivided leapfrog step.

ADAPTIVE FULL LEAPFROG STEP: r(n), v(n) –> r(n+1), v(n+1)

Physics-based time integration method with adaptive step refinement:

Subdivide the time step h = ΔT in powers-of-2 “micro-steps”
Compare a (2N+1)-step “coarse” vs. a (4N+1)-step “fine” integration
If within tolerance, return one of {coarse, fine, rich (Richardson extrapolation)}
Otherwise, double N and repeat until convergence or max subdivision reached

This helper function implements the leapfrog integration using a fixed number of micro-steps (subSteps, either 2N+1 for coarse or 4N+1 for fine). It takes the current radius and velocity (r_in, v_in) and applies a Kick-Drift-Kick sequence:

Initial half-kick using halfKick timestep.
A sequence of (subSteps - 1) / 2 - 1 full Drift-Kick pairs using midStep timestep.
Final full Drift using midStep.
Final half-kick using halfKick. The timestep sizes halfKick and midStep depend on whether it’s a coarse (h/(2N) and h/N) or fine (h/(4N) and h/(2N)) integration sequence. The function returns the new radius and velocity (r_out, v_out).

Parameters:

i – Particle index (used for rank in gravitational force).
npts – Total number of particles.
r_in – Input radius.
v_in – Input velocity.
ell – Angular momentum.
h – Total timestep size for the full adaptive step (ΔT).
N – Base subdivision factor for this micro-step sequence.
subSteps – Number of substeps (either 2N+1 or 4N+1).
grav – Gravitational constant value.
r_out – Pointer to store the output radius.
v_out – Pointer to store the output velocity.

static void doAdaptiveFullLeap(int i, int npts, double r_in, double v_in, double ell, double h, double radius_tol, double velocity_tol, int max_subdiv, double grav, int out_type, double *r_out, double *v_out)

Performs an adaptive full leapfrog step with error control.

This function implements the adaptive leapfrog algorithm to advance a particle’s state (r, v) over a full timestep h (ΔT). It iteratively refines the integration by comparing coarse (2N+1 substeps) and fine (4N+1 substeps) integrations performed via doMicroLeapfrog. The process starts with subdivision factor N=1. In each iteration:

A “coarse” pass using 2N+1 micro-steps is performed.
A “fine” pass using 4N+1 micro-steps is performed.
The relative differences between the final radius and velocity from the coarse and fine passes are compared against radius_tol and velocity_tol.
If both tolerances are met, the integration has converged. The final state (r_out, v_out) is determined based on out_type:
- 0: Use the coarse result.
- 1: Use the fine result.
- 2: Use Richardson extrapolation (4*fine - 3*coarse) for a potentially higher-order result.
If tolerances are not met, the subdivision factor N is doubled (N *= 2), and the process repeats, up to a maximum subdivision factor max_subdiv. If convergence is not achieved within max_subdiv, the function returns the result from the highest N iteration based on out_type.

See also

doMicroLeapfrog

Parameters:

i – Particle index (for rank in force computation).
npts – Total number of particles.
r_in – Input radius at the start of the step.
v_in – Input velocity at the start of the step.
ell – Angular momentum (conserved).
h – Full physical timestep size ΔT.
radius_tol – Relative convergence tolerance for radius.
velocity_tol – Relative convergence tolerance for velocity.
max_subdiv – Maximum subdivision factor N allowed.
grav – Gravitational constant value.
out_type – Result selection: 0=coarse, 1=fine, 2=Richardson extrapolation.
r_out – Pointer to store the final radius.
v_out – Pointer to store the final velocity.

static inline double dRhoDtaufun(double rhoVal, double vVal)

Calculates the derivative of rho with respect to the fictitious time tau.

LEVI-CIVITA REGULARIZATION

Physics-based regularized time integration method:

Transforms coordinates (r -> ρ = √r) to handle close encounters
Uses fictitious time τ to integrate equations of motion
Maps back to physical coordinates and time after integration
Provides enhanced stability for high-eccentricity orbits Part of the Levi-Civita regularization transformation.

Parameters:

rhoVal – The transformed radial coordinate ρ = √r
vVal – The velocity in transformed coordinates

Returns:

dρ/dτ = 0.5 * ρ * v

static inline double forceLCfun(int i, int npts, double totalmass, double grav, double ell, double rhoVal)

Calculates the total force acting on a particle in Levi-Civita transformed coordinates.

This function combines the gravitational and angular components of force in the ρ (rho) coordinate system which is the square root transformation of radius.

Parameters:

i – Particle index for force computation
npts – Total number of particles
totalmass – Total mass of the system including dark matter halo
grav – Gravitational constant
ell – Angular momentum
rhoVal – Current value of ρ = √r in Levi-Civita coordinates

Returns:

Total force in transformed coordinates

static void doLeviCivitaLeapfrog(int i, int npts, double r_in, double v_in, double ell, double dt, int N_taumin, double grav, double *r_out, double *v_out)

Performs integration of particle motion using Levi-Civita regularization.

This function implements the Levi-Civita leapfrog integration algorithm to advance a particle’s state over a physical time step dt.

Transforms coordinates from physical r_in to regularized rho = sqrt(r_in).
Estimates an initial fictitious time step deltaTau aiming for roughly N_taumin steps within the physical interval dt.
Integrates the equations of motion in (rho, v_rad, t_phys) using a leapfrog scheme in fictitious time tau:
- Compute force fval at rho_cur.
- Half-kick velocity: v_half = v_cur + 0.5 * deltaTau * fval.
- Drift rho: rho_next = rho_cur + deltaTau * dRhoDtaufun(rho_cur, v_half).
- Compute force fval2 at rho_next.
- Second half-kick velocity: v_next = v_half + 0.5 * deltaTau * fval2.
- Update physical time t_phys using the midpoint rho: t_next = t_phys + deltaTau * rho_mid^2.
The loop continues until t_phys >= dt.
If a step overshoots dt (t_next > dt), the final state (rho_f, v_f) at exactly t_phys = dt is found using linear interpolation between the state before and after the overshoot step.
Transforms the final regularized state (rho_f, v_f) back to physical coordinates r_out = rho_f^2 and v_out = v_f. Includes a maximum step count to prevent infinite loops.

See also

forceLCfun

See also

dRhoDtaufun

Parameters:

i – Particle index (for rank in force calculation).
npts – Total number of particles.
r_in – Initial physical radius at time t_i.
v_in – Initial physical radial velocity at time t_i.
ell – Angular momentum for the force law.
dt – Physical time step ΔT to advance.
N_taumin – Target number of tau-steps (influences deltaTau guess).
grav – Gravitational constant value.
r_out – Pointer to store the final physical radius at time t_i + ΔT.
v_out – Pointer to store the final physical velocity at time t_i + ΔT.

static void doMicroLeviCivita(int i, int npts, double rho_in, double v_in, double t_in, int subSteps, double h_tau, double grav, double ell, double *rho_out, double *v_out, double *t_out): Performs Levi-Civita regularized integration for a fixed number of micro-steps.

ADAPTIVE FULL LEVI-CIVITA REGULARIZATION

Enhanced regularization scheme that combines adaptive step sizing with Levi-Civita coordinate transformation for optimal performance near the coordinate origin.

Integrates the equations of motion in (rho, v, t) using fictitious time tau for a specified number of subSteps. This is a helper function used by doSingleTauStepAdaptiveLeviCivita.

Parameters

i : int Particle index (for rank in force calculation). npts : int Total number of particles in the simulation. rho_in : double Initial rho = sqrt(r) at the start of the tau interval h_tau. v_in : double Initial radial velocity at the start of the tau interval. t_in : double Initial physical time at the start of the tau interval. subSteps : int Number of micro-steps to perform (e.g., 2N+1 or 4N+1). h_tau : double Total fictitious time interval (Delta tau) for this integration sequence. grav : double Gravitational constant value (e.g., G_CONST). ell : double Angular momentum per unit mass (conserved). rho_out : double* Pointer to store the final rho value after h_tau. v_out : double* Pointer to store the final velocity after h_tau. t_out : double* Pointer to store the final physical time after h_tau.

Returns

None (results are returned via rho_out, v_out, t_out pointers).

See also

doSingleTauStepAdaptiveLeviCivita

See also

forceLCfun

See also

dRhoDtaufun

static void doSingleTauStepAdaptiveLeviCivita(int i, int npts, double rho_in, double v_in, double t_in, double h_guess, double radius_tol, double velocity_tol, int max_subdiv, double grav, double ell, int out_type, double *rho_out, double *v_out, double *t_out)

Performs a single adaptive step in Levi-Civita coordinates with error control.

Integrates over a proposed fictitious time step h_guess (Δτ) using adaptive refinement. It compares coarse (2N+1 substeps) and fine (4N+1 substeps) integrations via doMicroLeviCivita. The process starts with N=1. In each iteration:

Perform coarse integration (2N+1 substeps) to get (rhoC, vC, tC).
Perform fine integration (4N+1 substeps) to get (rhoF, vF, tF).
Compare relative differences in rho and v against tolerances.
If converged: Return result (rho_out, v_out, t_out) based on out_type.
- out_type=0: Coarse result (rhoC, vC, tC).
- out_type=1: Fine result (rhoF, vF, tF).
- out_type=2: Richardson extrapolation result (4*fine - 3*coarse) applied to rho, v, and t for a higher-order estimate. A floor of 0.0 is applied to the extrapolated rho.
If not converged: Double N (N *= 2) and repeat, up to max_subdiv. If convergence is not achieved, the result from the highest N fine step (4*N+1 substeps) is calculated and returned.

See also

doMicroLeviCivita

Parameters:

i – Particle index (for rank in force calculation).
npts – Total number of particles.
rho_in – Initial rho = sqrt(r) at the start of this adaptive tau step.
v_in – Initial radial velocity at the start of this step.
t_in – Initial physical time at the start of this step.
h_guess – Proposed fictitious time step (Δτ guess) for this adaptive step.
radius_tol – Relative convergence tolerance for rho comparison.
velocity_tol – Relative convergence tolerance for velocity comparison.
max_subdiv – Maximum allowed value for the subdivision factor N.
grav – Gravitational constant value (e.g., G_CONST).
ell – Angular momentum per unit mass (conserved).
out_type – Result selection mode: 0=coarse, 1=fine, 2=Richardson extrapolation.
rho_out – Pointer to store the final rho value after the adaptive step.
v_out – Pointer to store the final velocity after the adaptive step.
t_out – Pointer to store the final physical time after the adaptive step.

Returns:

None (results are returned via rho_out, v_out, t_out pointers).

static void doAdaptiveFullLeviCivita(int i, int npts, double r_in, double v_in, double ell, double dt, int N_taumin, double radius_tol, double velocity_tol, int max_subdiv, double grav, int out_type, double *r_out, double *v_out)

Performs an adaptive full Levi-Civita regularized leapfrog step over physical time dt.

Integrates particle motion over a physical time step dt (ΔT) using the adaptive Levi-Civita scheme. It takes multiple adaptive tau steps (using doSingleTauStepAdaptiveLeviCivita) until the accumulated physical time t_cur reaches or exceeds dt.

Handles the edge case of near-zero initial radius r_in by returning immediately.
Transforms initial physical state (r_in, v_in) to regularized (rho_current, v_current).
Estimates an initial fictitious time step deltaTau based on r_in and N_taumin.
Enters a loop that continues as long as t_cur < dt: a. Calls doSingleTauStepAdaptiveLeviCivita to perform one adaptive τ-step, advancing the state to (rho_next, v_next, t_next). b. Checks if the step overshot the target time (t_next > dt). c. If overshoot: Performs linear interpolation between the state before (rho_current, v_current, t_cur) and after (rho_next, v_next, t_next) the step to find the state (rho_final, v_final) exactly at t = dt. Transforms back to physical (r_out, v_out) and returns. d. If no overshoot: Updates the current state (rho_current = rho_next, etc.) and continues the loop.
Includes a maximum step count (stepMax) as a safety measure against infinite loops.
If the loop completes normally (e.g., t_cur exactly equals dt on step acceptance, or if the overshoot check somehow fails), it transforms the final state (rho_current, v_current) back to physical coordinates (r_out, v_out).

See also

doSingleTauStepAdaptiveLeviCivita

Parameters:

i – Particle index (for rank in force calculation).
npts – Total number of particles.
r_in – Initial physical radius at the start of the physical step dt.
v_in – Initial physical radial velocity at the start of the step dt.
ell – Angular momentum per unit mass (conserved).
dt – Physical time step (ΔT) to advance.
N_taumin – Target number of tau-steps (influences deltaTau guess).
radius_tol – Relative convergence tolerance for rho comparison (passed to substep function).
velocity_tol – Relative convergence tolerance for velocity comparison (passed to substep function).
max_subdiv – Maximum allowed value for the subdivision factor N (passed to substep function).
grav – Gravitational constant value (e.g., G_CONST).
out_type – Result selection mode for substeps: 0=coarse, 1=fine, 2=Richardson.
r_out – Pointer to store the final physical radius after time dt.
v_out – Pointer to store the final physical velocity after time dt.

Returns:

None (results are returned via r_out, v_out pointers).

static void write_initial_conditions(double **particles, int npts, const char *filename): Writes particle initial conditions to a binary file.

INITIAL CONDITION FILE I/O FUNCTIONS

Stores the complete initial particle state (radius, velocity, angular momentum, original index, orientation parameter mu) and the particle count (npts) to a binary file for later retrieval via read_initial_conditions. Opens the file in write binary mode (“wb”). First writes the integer npts, then writes the 5 double-precision values for each particle sequentially.

Parameters

particles : double** 2D array containing particle properties [component][particle_index]. Expected components: 0=rad, 1=vel, 2=angmom, 3=orig_idx, 4=mu. npts : int Number of particles to write. filename : const char* Path to the output binary file.

Returns

None

See also

read_initial_conditions

Note

The binary file format is: npts (int32), followed by npts records, each consisting of 5 double values (radius, velocity, ang. mom., orig. index, mu).

Warning

Prints an error to stderr if the file cannot be opened.

static void read_initial_conditions(double **particles, int npts, const char *filename)

Reads particle initial conditions from a binary file.

Loads the complete particle state (radius, velocity, angular momentum, original index, orientation parameter mu) from a binary file previously created by write_initial_conditions. Verifies that the number of particles read from the file matches the expected count npts. Opens the file in read binary mode (“rb”).

const char *get_sort_description(const char *sort_alg)

Returns a human-readable description for a sort algorithm identifier string.

Maps internal sort algorithm identifiers to user-friendly descriptive names. Supported algorithms: quadsort_parallel, quadsort, insertion_parallel, insertion.

Parameters:: sort_alg – The internal algorithm identifier string (e.g., “quadsort_parallel”)
Returns:: A descriptive name for the algorithm (e.g., “Parallel Quadsort”), or the input string itself if no match is found

static void parseSaveArgs(int argc, char *argv[], int *pIndex)

Parses sub-arguments for the “–save” command-line option.

Reads subsequent arguments after “–save” (until another option starting with ‘-’ is encountered) and sets global data output flags (g_doDebug, g_doDynPsi, g_doDynRank, g_doAllParticleData) based on the highest priority sub-argument found. Valid sub-arguments: “raw-data”, “psi-snaps”, “full-snaps”, “debug-energy”, “all”.

static void reassign_orig_ids_with_rank(double *orig_ids, int n)

Remaps original particle IDs to their rank (0…n-1) among the final set.

Takes an array of potentially non-contiguous original particle IDs (stored as doubles) and replaces each ID with its zero-based rank within the sorted sequence of those IDs. This effectively transforms arbitrary ID values into a contiguous sequence [0, n-1]. Used after tidal stripping.

void fft_gaussian_convolution(const double *density_grid, int grid_size, const double *log_r_grid, double sigma_log, double *result)

Applies Gaussian smoothing using FFT-based convolution (thread-safe via critical section).

SIGNAL PROCESSING AND FILTERING UTILITIES

Advanced numerical processing utilities for density field handling including:

FFT-based convolution for density smoothing
Direct Gaussian convolution for smaller datasets
Signal filtering and processing functions

FFT METHODS AND CONVOLUTION IMPLEMENTATIONS

Smooths a density field defined on a potentially non-uniform grid (log_r_grid) using FFT convolution with a Gaussian kernel of width sigma_log (defined in log-space). Uses zero-padding to avoid wrap-around artifacts. This is generally faster than direct convolution for large grid_size. Assumes log_r_grid is uniformly spaced.

Parameters

density_grid : const double* Input density grid array (values corresponding to log_r_grid). grid_size : int Number of points in the input grid and density arrays. log_r_grid : const double* Array of logarithmic radial grid coordinates (log10(r)). Must be uniformly spaced. sigma_log : double Width (standard deviation) of the Gaussian kernel in log10-space. result : double* Output array (pre-allocated, size grid_size) for the smoothed density field.

Returns

None (populates the result array).

See also

direct_gaussian_convolution

See also

gaussian_convolution

Note

Uses FFTW library for Fast Fourier Transforms (fftw_malloc, fftw_plan_dft_r2c_1d, etc.). Requires FFTW to be installed. FFTW operations are protected by omp critical(fftw).

Note

Resulting smoothed density is clamped to a minimum value of 1e-10.

Warning

Prints errors to stderr and returns early on memory allocation failures or FFTW plan creation failures.

void direct_gaussian_convolution(const double *density_grid, int grid_size, const double *log_r_grid, double sigma_log, double *result)

Applies Gaussian smoothing using direct convolution.

Smooths a density field defined on a potentially non-uniform grid (log_r_grid) using direct spatial convolution with a Gaussian kernel of width sigma_log (defined in log-space). For each point i in the output result array, it computes a weighted sum of the input density_grid values: result[i] = sum(density_grid[j] * kernel(log_r_grid[i] - log_r_grid[j])) / sum(kernel(...)) where the kernel is a Gaussian function G(x) = (1/(σ√2π)) * exp(-0.5*(x/σ)²), with x being the distance log_r_grid[i] - log_r_grid[j] and σ = sigma_log. This method is generally more accurate than FFT-based convolution, especially for non-uniform grids or near boundaries, but has a higher computational cost (O(N²)) which makes it slower for large grid_size.

void gaussian_convolution(const double *density_grid, int grid_size, const double *log_r_grid, double sigma_log, double *result)

Performs Gaussian smoothing by selecting the appropriate convolution method.

Acts as a routing function that delegates density smoothing to either direct_gaussian_convolution or fft_gaussian_convolution based on the value of the global debug_direct_convolution flag.

If debug_direct_convolution is non-zero, direct_gaussian_convolution is called (more accurate, O(N²) complexity, suitable for smaller or non-uniform grids).
If debug_direct_convolution is zero (default), fft_gaussian_convolution is called (faster for large grids, O(N log N) complexity, requires uniformly spaced logarithmic grid).

static int find_last_processed_snapshot(int *snapshot_steps, int noutsnaps)

Finds the index of the last successfully processed snapshot in restart mode.

RESTART AND RECOVERY MANAGEMENT

Checks for the existence and basic validity (non-zero size, consistent size compared to snapshot 0) of expected output data files (unsorted and sorted Rank/Mass/Rad/Vrad files) for each snapshot step defined in snapshot_steps. Used to determine where to resume processing when the --restart flag is enabled.

Parameters

snapshot_steps : int* Array containing the step numbers (indices) of the snapshots that were supposed to be written during the simulation. noutsnaps : int Total number of snapshot indices in the snapshot_steps array.

Returns

int

The index s in snapshot_steps corresponding to the next snapshot that needs processing (i.e., last_valid_index + 1).
Returns -1 if no valid snapshot files are found (start from beginning).
Returns -2 if all expected snapshot files exist and appear valid (nothing to do).

Note

Compares file sizes against the first snapshot’s files as a reference, allowing for a +/- 5% tolerance. Assumes g_file_suffix is set correctly.

Warning

Prints status messages and warnings/errors to stdout/log.

int main(int argc, char *argv[])

Main entry point for the n-sphere dark matter simulation program.

Orchestrates the overall simulation workflow:

Parses command-line arguments using a flexible option system
Sets up initial conditions (particle distribution, system parameters)
Executes the selected integration method for trajectory evolution
Outputs data products (particle states, phase diagrams, energy tracking)
Handles restart/resume functionality when requested

double massintegrand (double r, void *params __attribute__((unused)))

Calculates the mass integrand for a given radius using a density profile.

Parameters:

r – Radius value
params – Unused parameter (with compiler attribute to prevent warnings)

Returns:

Mass integrand value at the specified radius

static void insertion_sort(double **columns, int n)

Implementation of classic insertion sort algorithm for particle data.

Sorts an array of particle data using the insertion sort algorithm. While not the fastest algorithm for large datasets, it is stable and works well for small arrays or nearly sorted data.

Parameters:

columns – 2D array of particle data to be sorted
n – Number of elements to sort

static void stdlib_qsort_wrapper(double **columns, int n)

Wrapper for the standard C library quicksort function.

Provides a consistent interface to the standard library qsort function using the compare_particles function as the comparison callback.

Parameters:

columns – 2D array of particle data to be sorted (passed as double**).
n – Number of elements (particles) to sort.

static void quadsort_wrapper(double **columns, int n)

Wrapper for the external quadsort algorithm.

Provides a consistent interface to the external quadsort function, which is typically faster than standard quicksort for many distributions.

Parameters:

columns – 2D array of particle data to be sorted
n – Number of elements to sort

static void insertion_sort_sub(double **columns, int start, int end)

Helper function that performs insertion sort on a subarray of particle data.

Similar to insertion_sort but operates on a specific range [start..end] inclusive. Used by the parallel sorting algorithms to sort individual chunks and seam regions.

Parameters:

columns – 2D array of particle data containing the target subarray (double**).
start – Starting index of the subarray (inclusive).
end – Ending index of the subarray (inclusive).

void insertion_parallel_sort(double **columns, int n)

Parallel insertion sort implementation using chunk-based approach with overlap.

Divides the data into num_sort_sections chunks, sorts each chunk in parallel, then fixes the boundaries between chunks by re-sorting overlap regions. This approach balances parallelism with the need to ensure properly sorted output.

Parameters:

columns – 2D array of particle data to be sorted
n – Number of elements to sort

void quadsort_parallel_sort(double **columns, int n)

Parallel quadsort implementation using chunk-based approach with overlap.

Similar to insertion_parallel_sort but uses the quadsort algorithm for both the initial chunk sorting and the seam fixing. Quadsort is generally faster than insertion sort for larger datasets while maintaining stability.

Parameters:

columns – 2D array of particle data to be sorted
n – Number of elements to sort

void verify_sort_results(double **columns, int n, const char *label)

Validates sorting results by comparing with standard qsort.

Creates a copy of the input array, sorts it with standard qsort, then compares the results element by element to verify that the sorting algorithm produced the expected ordering.

Parameters:

columns – 2D array of particle data that has been sorted
n – Number of elements in the array
label – Name of the sorting algorithm for diagnostic output

int compare_partdata_by_rad(const void *a, const void *b)

Comparison function for sorting PartData structures by radius value.

Provides a robust comparison that handles NULL pointers and NaN values, ensuring stable sorting behavior even with potentially problematic data.

Parameters:

a – Pointer to the first PartData structure (as void pointer)
b – Pointer to the second PartData structure (as void pointer)

Returns:

-1 if a<b, 1 if a>b, 0 if equal

int compare_by_rr(const void *a, const void *b)

Comparison function for sorting RrPsiPair structures by radius (rr).

See also

RrPsiPair

See also

sort_rr_psi_arrays

Parameters:

a – Pointer to the first RrPsiPair structure (as void pointer).
b – Pointer to the second RrPsiPair structure (as void pointer).

Returns:

-1 if a->rr < b->rr, 1 if a->rr > b->rr, 0 if equal.

Variables

static char g_file_suffix[256] = "": Global file suffix string.

int g_enable_logging = 0

Writes a formatted message to the log file with timestamp and severity level.

Enable logging to file (controlled by —enable-log flag).

int g_doDebug = 0

Global simulation feature flags.

These control various optional behaviors and optimizations in the simulation. Each flag can be set via command line arguments. Enable detailed debug output.

int g_doDynPsi = 0: Enable dynamic potential recalculation.

int g_doDynRank = 0: Enable dynamic rank calculation per step.

int g_doAllParticleData = 0: Save complete particle evolution history.

int g_doRestart = 0: Enable simulation restart from checkpoint.

int skip_file_writes = 0: Skip file writes during simulation restart.

int debug_direct_convolution = 0

Convolution method selection for density smoothing.

Controls the convolution method selection for density smoothing.

0 = FFT-based convolution (default, faster for large datasets). 1 = Direct spatial convolution (more accurate but slower).

Global flag that determines whether to use direct convolution (1) or FFT-based convolution (0). Direct convolution is more accurate but slower for large grids, while FFT-based convolution is faster but may introduce artifacts.

See also

gaussian_convolution

Note

External declaration, defined elsewhere.

static int use_closest_to_Lcompare = 1

Angular momentum selection configuration for particle filtering.

Mode 0: Select particles with the 5 lowest L values. Mode 1: Select particles with L values closest to Lcompare. Mode selector (0 or 1).

static double Lcompare = 0.05: Reference L value for closest-match mode (Mode 1).

static int *Rank_partdata_snap = NULL: Global particle data arrays for snapshot processing.

GLOBAL PARTICLE DATA ARRAYS

Particle rank (sorted position) data for a snapshot.

static float *R_partdata_snap = NULL: Radial position data for a snapshot.

static float *Vrad_partdata_snap = NULL: Radial velocity data for a snapshot.

static float *L_partdata_snap = NULL: Angular momentum data for a snapshot.

static float *L_block = NULL

Global arrays for particle data processing (block storage).

Angular momentum block.

static int *Rank_block = NULL: Particle rank (sorted position) block.

static float *R_block = NULL: Radial position block.

static float *Vrad_block = NULL: Radial velocity block.

static int nlowest = 5

Variables for tracking low angular momentum particles.

Number of lowest angular momentum particles to track.

static int *chosen = NULL: Array of indices (original IDs) for selected low-L particles.

static double **lowestL_r = NULL: Radial positions of tracked low-L particles over time [particle][time_step].

static double **lowestL_E = NULL: Energy values of tracked low-L particles over time [particle][time_step].

static double **lowestL_L = NULL: Angular momenta of tracked low-L particles over time [particle][time_step].

static int use_identity_gravity = 0

Gravitational force calculation control flag.

Used for testing and debugging orbital dynamics:

0 = Normal gravitational force calculation (default).
1 = Zero gravity (particles move in straight lines).

Note

Setting this to 1 is useful for validating the integration scheme independent of gravitational physics.

static double dbg_approxE[DEBUG_MAX_STEPS]

Arrays for tracking energy components through simulation for debugging.

Theoretical model energy (per unit mass).

static double dbg_dynE[DEBUG_MAX_STEPS]: Actual dynamical energy (per unit mass).

static double dbg_kinE[DEBUG_MAX_STEPS]: Kinetic energy component (per unit mass).

static double dbg_potE[DEBUG_MAX_STEPS]: Potential energy component (per unit mass).

static double dbg_time[DEBUG_MAX_STEPS]: Simulation time at each snapshot (Myr).

static double dbg_radius[DEBUG_MAX_STEPS]: Particle radius at each snapshot (kpc).

static int dbg_count = 0: Number of debug snapshots recorded.

static double normalization: Mass normalization factor for energy calculations. Calculated based on the integral of the density profile.

static int doReadInit = 0

Flag indicating whether to read initial conditions from file (1=yes, 0=no).

FILE I/O AND DATA MANAGEMENT SUBSYSTEM

Functions for saving, loading, and managing simulation data including:

Initial condition generation and I/O
Snapshot file management
Binary file format utilities

static int doWriteInit = 0: Flag indicating whether to write initial conditions to file (1=yes, 0=no).

static const char *readInitFilename = NULL: Filename to read initial conditions from (if doReadInit=1).

static const char *writeInitFilename = NULL: Filename to write initial conditions to (if doWriteInit=1).

static const char *g_defaultSortAlg = "quadsort_parallel": < Default sorting algorithm identifier string. Set based on command-line options.

SORTING ALGORITHM CONFIGURATION

static const int num_sort_sections = 24: Constants and functions controlling the parallel sorting algorithm behavior.

PARALLEL SORTING ALGORITHM CONSTANTS AND FUNCTIONS

The sorting implementation uses a parallel chunk-based approach with overlapping regions between chunks to ensure correct ordering at chunk boundaries. Number of parallel chunks for sorting.

static const int overlap_divisor = 8: Overlap factor: Overlap = section_size / overlap_divisor.

static const int min_seam_overlap = 50: Minimum overlap region size between chunks.

struct fEintegrand_params

Parameters for energy integration calculations.

ENERGY CALCULATION AND INTEGRATION STRUCTURES

Used as void* params argument in GSL integration routines, specifically for distribution function calculations (fEintegrand).

Public Members

double E: Energy value (relative energy).

gsl_spline *splinePsi: Interpolation spline for potential Psi(r).

gsl_spline *splinemass: Interpolation spline for enclosed mass M(r).

gsl_interp_accel *rofPsiarray: Accelerator for radius lookup from potential r(Psi).

gsl_interp_accel *massarray: Accelerator for mass lookups M(r).

struct LAndIndex

Structure to track angular momentum with particle index and direction.

Used in sorting and selection of particles by angular momentum, especially when finding particles closest to a reference L.

Public Members

double L: Angular momentum value (or squared difference from Lcompare).

int idx: Original particle index (before sorting by L).

int sign: Direction indicator (+1 or -1) or sign of (L - Lcompare).

struct PartData

Compact data structure for particle properties used in sorting and analysis.

PARTICLE DATA STRUCTURES AND OPERATIONS

Contains essential physical properties (radial position, velocity, angular momentum) and tracking metadata (rank, original index) for each particle in the simulation. Used extensively for sorting, file I/O, and data analysis operations.

Note

Uses compact float types for physical quantities to reduce memory usage when processing large particle counts.

Note

The rank field is assigned during radial sorting, while original_index preserves the initial array position for tracking particles across snapshots.

Public Members

int rank

float rad

float vrad

float angmom

int original_index

struct RrPsiPair

Utility functions and structures for sorting spline data arrays.

SPLINE DATA SORTING UTILITIES

Provides mechanisms to sort arrays used for GSL spline creation (like radius r and potential Psi) while maintaining the correct correspondence between paired values after sorting based on one of the arrays (typically radius).

Structure for pairing radius and potential values during sorting.

Used as a temporary structure within sort_rr_psi_arrays to maintain the correspondence between radius and potential values when sorting by radius.

Public Members

double rr: Radius value.

double psi: Corresponding potential value.

NSphere Sort (nsphere_sort.h)

Clean wrapper for quadsort to suppress warnings and add support for doubles.

This header provides a clean interface to quadsort by:

Containing all warnings through pragma directives
Adding support for double and long double types
Keeping the original quadsort implementation pristine

Note

This wrapper is specifically designed for the NSphere simulation system to ensure consistent sorting behavior across all numerical types.

Defines

ns_sort_array

Restore warning settings to their previous state.

Pops the diagnostic pragma stack to restore compiler warning settings that were modified earlier in this header.

Define clean aliases to quadsort functions with zero overhead

These macros provide a consistent naming convention for sorting functions in the NSphere codebase, making it clear which functions are part of the NSphere API rather than direct quadsort calls.

ns_sort_prim

ns_sort_size

ns_sort_extended

Functions

static inline void ns_quadsort_extended(void *array, size_t nmemb, size_t size, CMPFUNC *cmp)

Push current warning state and disable specific warnings.

These pragma directives are used to suppress compiler warnings that would otherwise be generated by the quadsort implementation. This allows us to use the original quadsort code without modification.

NSphere extensions to quadsort - Support for double and long double types

These extensions add specialized sorting implementations for floating-point types that may not be efficiently handled by the standard quadsort. They are conditionally compiled based on compiler support.

Extended quadsort function with support for double and long double types

This function extends the standard quadsort implementation to handle double and long double types with specialized implementations. For other types, it falls back to the standard quadsort. This provides optimal performance for all data types used in the NSphere simulation system.

QuadSort (quadsort.c / quadsort.h)

Functions

void FUNC() parity_swap_four (VAR *array, CMPFUNC *cmp)

void FUNC() parity_swap_five (VAR *array, CMPFUNC *cmp)

void FUNC() parity_swap_six (VAR *array, VAR *swap, CMPFUNC *cmp)

void FUNC() parity_swap_seven (VAR *array, VAR *swap, CMPFUNC *cmp)

void FUNC() tiny_sort (VAR *array, VAR *swap, size_t nmemb, CMPFUNC *cmp)

void FUNC() parity_merge (VAR *dest, VAR *from, size_t left, size_t right, CMPFUNC *cmp)

void FUNC() tail_swap (VAR *array, VAR *swap, size_t nmemb, CMPFUNC *cmp)

void FUNC() quad_reversal (VAR *pta, VAR *ptz)

void FUNC() quad_swap_merge (VAR *array, VAR *swap, CMPFUNC *cmp)

void FUNC() tail_merge (VAR *array, VAR *swap, size_t swap_size, size_t nmemb, size_t block, CMPFUNC *cmp)

size_t FUNC() quad_swap (VAR *array, size_t nmemb, CMPFUNC *cmp)

void FUNC() cross_merge (VAR *dest, VAR *from, size_t left, size_t right, CMPFUNC *cmp)

void FUNC() quad_merge_block (VAR *array, VAR *swap, size_t block, CMPFUNC *cmp)

size_t FUNC() quad_merge (VAR *array, VAR *swap, size_t swap_size, size_t nmemb, size_t block, CMPFUNC *cmp)

void FUNC() partial_forward_merge (VAR *array, VAR *swap, size_t swap_size, size_t nmemb, size_t block, CMPFUNC *cmp)

void FUNC() partial_backward_merge (VAR *array, VAR *swap, size_t swap_size, size_t nmemb, size_t block, CMPFUNC *cmp)

void FUNC() trinity_rotation (VAR *array, VAR *swap, size_t swap_size, size_t nmemb, size_t left)

size_t FUNC() monobound_binary_first (VAR *array, VAR *value, size_t top, CMPFUNC *cmp)

void FUNC() rotate_merge_block (VAR *array, VAR *swap, size_t swap_size, size_t lblock, size_t right, CMPFUNC *cmp)

void FUNC() rotate_merge (VAR *array, VAR *swap, size_t swap_size, size_t nmemb, size_t block, CMPFUNC *cmp)

void FUNC() quadsort (void *array, size_t nmemb, CMPFUNC *cmp)

void FUNC() quadsort_swap (void *array, void *swap, size_t swap_size, size_t nmemb, CMPFUNC *cmp)

Defines

QUAD_CACHE

head_branchless_merge(ptd, x, ptl, ptr, cmp)

tail_branchless_merge(tpd, y, tpl, tpr, cmp)

parity_merge_two(array, swap, x, ptl, ptr, pts, cmp)

parity_merge_four(array, swap, x, ptl, ptr, pts, cmp)

branchless_swap(pta, swap, x, cmp)

swap_branchless(pta, swap, x, y, cmp)

VAR

FUNC(NAME)

VAR

FUNC(NAME)

cmp(a, b)

VAR

FUNC(NAME)

cmp(a, b)

VAR

FUNC(NAME)

VAR

FUNC(NAME)

cmp(a, b)

VAR

FUNC(NAME)

cmp(a, b)

QUAD_CACHE

VAR

FUNC(NAME)

VAR

FUNC(NAME)

Typedefs

typedef int CMPFUNC(const void *a, const void *b)

Functions

void quadsort(void *array, size_t nmemb, size_t size, CMPFUNC *cmp)

void quadsort_prim(void *array, size_t nmemb, size_t size)

void quadsort_size(void *array, size_t nmemb, size_t size, CMPFUNC *cmp)

C API Reference

NSphere Core (nsphere.c)

GLOBAL CONFIGURATION AND MACROS

PHYSICAL CONSTANTS AND ASTROPHYSICAL PARAMETERS

Energy Debugging and Validation Subsystem

BINARY FILE I/O UTILITIES

Parameters

Returns

PARTICLE SORTING AND ORDERING OPERATIONS

Parameters

Returns

PHYSICS CALCULATION FUNCTIONS

Parameters

Returns

COMMAND LINE ARGUMENT PROCESSING

Parameters

Returns

ADAPTIVE FULL LEAPFROG STEP: r(n), v(n) –> r(n+1), v(n+1)

LEVI-CIVITA REGULARIZATION

ADAPTIVE FULL LEVI-CIVITA REGULARIZATION

Parameters

Returns

INITIAL CONDITION FILE I/O FUNCTIONS

Parameters

Returns

SIGNAL PROCESSING AND FILTERING UTILITIES

FFT METHODS AND CONVOLUTION IMPLEMENTATIONS

Parameters

Returns

RESTART AND RECOVERY MANAGEMENT

Parameters

Returns

GLOBAL PARTICLE DATA ARRAYS

FILE I/O AND DATA MANAGEMENT SUBSYSTEM

SORTING ALGORITHM CONFIGURATION

PARALLEL SORTING ALGORITHM CONSTANTS AND FUNCTIONS

ENERGY CALCULATION AND INTEGRATION STRUCTURES

PARTICLE DATA STRUCTURES AND OPERATIONS

SPLINE DATA SORTING UTILITIES

NSphere Sort (nsphere_sort.h)

QuadSort (quadsort.c / quadsort.h)