NSphere Core (nsphere.c)

Complete API for the main NSphere simulation engine.

Defines

CLEAN_EXIT(code)

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

This macro ensures proper resource cleanup before program termination:

• Uses an OpenMP critical section to ensure only one thread performs cleanup.

• Calls cleanup_all_particle_data() 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.

CLEAN_LOCAL_EXIT(code)

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

Performs resource cleanup and exits in single-threaded contexts. This macro calls cleanup_all_particle_data() before exiting with the specified code.

Example

if (error_condition) {
    CLEAN_LOCAL_EXIT(1);
}

Warning

This macro is not thread-safe. Use it only in single-threaded sections of the code. For thread-safe exits, use CLEAN_EXIT.

Parameters:

• code -- The exit code for the program.

cube(x)

Calculates the cube of a value: \(x^3\).

CUTOFF_FACTOR_CORED_DEFAULT

Default \(r_{\text{max}}\) factor for Cored profile ( \(r_{\text{max}} = \text{factor} \times r_c\)).

CUTOFF_FACTOR_NFW_DEFAULT

Default \(r_{\text{max}}\) factor for NFW profile ( \(r_{\text{max}} = \text{factor} \times r_c\)).

DEBUG_MAX_STEPS

Maximum number of debug energy snapshots.

DEBUG_PARTICLE_ID

Advances azimuthal angle via 4th-order Simpson's rule using radius history.

Integrates \(d\phi/dt = L/r^2\) over one timestep using a three-point radius stencil at \(t-\Delta t\), \(t\), and \(t+\Delta t\). For the first timestep ( \(j=0\)), falls back to 2nd-order trapezoidal rule. For subsequent steps, uses quadratic interpolation to estimate \(r(t+\Delta t/2)\) and applies Simpson's rule: \(\Delta\phi = (\omega_t + 4\omega_{mid} + \omega_{t+\Delta t})\Delta t/6\). Updates \((\cos\phi, \sin\phi)\) via rotation matrix for exact norm preservation.

See also

effective_angular_force Particle ID (by rank) to track for energy debugging.

Note

Skips update if any radius \(\leq 10^{-15}\) or \(|L| \leq 10^{-15}\). Renormalizes \((\cos\phi, \sin\phi)\) if \(|\cos^2\phi + \sin^2\phi - 1| > 10^{-12}\).

Parameters:

• particles -- [in,out] Particle data array. Reads angular momentum [2][i] and updated position [0][i]. Updates azimuthal components [5][i] and [6][i].

• i -- [in] Particle index (0 to npts-1).

• dt -- [in] Physical timestep size (Myr).

• j -- [in] Current timestep index (0 triggers trapezoidal bootstrap).

• rlast -- [in] Radius at previous timestep \(r(t-\Delta t)\) (kpc).

• rcurrent -- [in] Radius at current timestep \(r(t)\) (kpc).

FALLOFF_FACTOR_NFW_DEFAULT

Default falloff concentration parameter \(C\) for NFW profile.

G_CONST

Newton's gravitational constant in kpc (km/sec) \(^2\)/ \(M_{\odot}\).

HALO_MASS

Default general halo mass (used for Cored profile by default) in \(M_{\odot}\).

HALO_MASS_NFW

Default NFW profile halo mass in solar masses ( \(M_{\odot}\)).

HIST_NBINS

Number of bins for 2D phase-space histograms (400x400 resolution).

IF_ALL_PART

Macro for code blocks executed only if g_doAllParticleData is true.

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.

imin(a, b)

Minimum of two integer values.

kmsec_to_kpcmyr

Conversion factor: km/s to kpc/Myr.

MYBINIO_H

NUM_MINI_SUBSTEPS_BOOTSTRAP

Number of mini Euler steps per bootstrap full step.

PI

Mathematical constant \(\pi\).

RADIX_BITS

Number of bits processed per radix sort pass (8 bits = 256 buckets).

RADIX_MASK

Bit mask for extracting radix bucket index (0xFF).

RADIX_SIZE

Number of buckets for radix sort ( \(2^8 = 256\)).

RC

Core radius in kpc.

RC_NFW_DEFAULT

Default NFW profile scale radius (kpc).

sqr(x)

Calculates the square of a value: \(x^2\).

VEL_CONV_SQ

Velocity conversion squared (kpc/Myr) \(^2\) per (km/s) \(^2\).

Typedefs

typedef double (*drhodr_func_t)(double r)

Function pointer types for dynamic dispatch.

These types enable runtime selection of profile-specific implementations without conditional branching in performance-critical loops. Function pointer for density derivative: \(d\rho/dr\) (Cored/Hernquist profiles).

typedef double (*drhodr_nfw_func_t)(double r, double rc, double nt_nfw, double falloff_C)

Function pointer for NFW density derivative: \(d\rho/dr\) with profile parameters.

typedef double (*potential_function_t)(double r, void *params)

Function pointer for potential evaluation: \(\Psi(r)\).

Functions

static int __attribute__ ((unused))

static void add_snapshot_to_double_buffer(double **particles, double *phi_array, int *rank_array, int npts)

Adds the current snapshot to the double-precision circular buffer.

Copies particle state (R, Vrad, L, Rank, phi) from current simulation state into the circular buffer at the current slot position. Updates buffer slot position and count. Buffer grows from 0 to 4 snapshots, then maintains rolling last 4.

Parameters:

• particles -- Pointer to particle state arrays [R, Vrad, L, ...]

• phi_array -- Array of particle phi angles

• rank_array -- Array of particle ranks (sorted positions)

• npts -- Number of particles

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

Adjusts the total number of timesteps to align with desired output snapshot intervals.

This function calculates an adjusted number of total simulation timesteps, \(N'_{times}\), such that it is greater than or equal to the initially requested Ntimes_initial ( \(N\)) and satisfies the constraint: \((N'_{times} - 1)\) must be an integer multiple of \((M - 1) \times p\). Here, \(M\) is nout (number of desired output snapshot points, which means \(M-1\) intervals) and \(p\) is dtwrite (the low-level write interval in terms of simulation timesteps). This alignment ensures that exactly nout snapshots can be produced at intervals that are multiples of dtwrite and that also evenly span the total adjusted simulation duration.

Parameters:

• Ntimes_initial -- [in] Initially requested total number of simulation timesteps ( \(N\)).

• nout -- [in] Number of desired output snapshot points ( \(M\)). Must be >= 2 for adjustment to apply.

• dtwrite -- [in] The interval (in timesteps) at which low-level data is potentially written ( \(p\)). Must be >= 1.

Returns:

int The adjusted total number of timesteps ( \(N'_{times}\)). Returns Ntimes_initial if nout < 2 or dtwrite < 1 or other edge cases where the constraint cannot be met.

static void allocate_double_precision_buffers(int npts)

Allocates memory for the double-precision snapshot buffer system.

Allocates circular buffers to store the last 4 snapshots in full double precision. This prevents precision loss during restart/extend operations by avoiding the double→float32→double conversion inherent in the regular all_particle_data files.

Buffer grows incrementally: stores [0], then [0,1], then [0,1,2], then [0,1,2,3], then rolls to maintain the last 4 snapshots.

See also

cleanup_all_particle_data()

Note

Exits via CLEAN_EXIT(1) if allocation fails.

Parameters:

npts -- Number of particles in the simulation.

static void allocate_trajectory_buffers(int num_traj_particles, int nlowest, int buffer_size, int npts)

Allocates memory for all trajectory buffering arrays.

This function is called once at the start of the simulation to allocate fixed-size buffers for storing trajectory data. The buffer size is determined by the snapshot writing schedule to ensure data is flushed to disk in synchronized blocks.

Parameters:

• num_traj_particles -- The number of low-ID particles to track.

• nlowest -- The number of low-L particles to track.

• buffer_size -- The total number of timesteps the buffer can hold.

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 full particle data to the specified output file.

This function writes a chunk of particle data, corresponding to block_size timesteps, to the given binary file. The data for each particle (rank, radius, radial velocity, angular momentum) is written sequentially for each timestep within the block (step-major order). This means all particle data for step s is contiguous, followed by all data for step s+1, etc., within the block. The file is opened in append binary mode ("ab"). This is used for creating the all_particle_data.dat file which stores the complete evolution history when g_doAllParticleData is enabled.

Note

Exits via CLEAN_EXIT(1) if the file cannot be opened for appending.

Parameters:

• filename -- [in] Path to the output binary file (e.g., "data/all_particle_data<suffix>.dat").

• npts -- [in] Number of particles.

• block_size -- [in] Number of timesteps of data contained in the provided _block arrays.

• L_block -- [in] Pointer to the block of angular momentum data (float array). Assumed to be [step_in_block * npts + particle_orig_id].

• Rank_block -- [in] Pointer to the block of particle rank data (int array). Assumed to be [step_in_block * npts + particle_orig_id].

• R_block -- [in] Pointer to the block of radial position data (float array). Assumed to be [step_in_block * npts + particle_orig_id].

• Vrad_block -- [in] Pointer to the block of radial velocity data (float array). Assumed to be [step_in_block * npts + particle_orig_id].

static void append_all_particle_ids_to_file(const char *filename, int npts, int block_size, int *id_block)

Appends a block of particle IDs to the specified output file.

This function writes a chunk of particle ID data, corresponding to block_size timesteps, to the given binary file. The data is written in step-major order to match the structure of all_particle_data.dat.

Parameters:

• filename -- [in] Path to the output binary ID data file.

• npts -- [in] Number of particles.

• block_size -- [in] Number of timesteps of data contained in id_block.

• id_block -- [in] Pointer to the block of particle ID data (int array).

static void append_all_particle_phi_data_chunk_to_file(const char *filename, int npts, int block_size, float *phi_block)

Appends a block of phi angle data to the specified output file.

This function writes a chunk of phi data, corresponding to block_size timesteps, to the given binary file. The data is written sequentially for each timestep within the block (step-major order). The file is opened in append binary mode ("ab"). This is used for creating the corresponding all_particle_data_phi.dat file.

Parameters:

• filename -- [in] Path to the output binary phi data file.

• npts -- [in] Number of particles.

• block_size -- [in] Number of timesteps of data contained in phi_block.

• phi_block -- [in] Pointer to the block of phi data (float array). Assumed to be [step_in_block * npts + particle_orig_id].

static void append_all_particle_scatter_counts_to_file(const char *filename, int npts, int block_size, int *scat_count_block)

Appends a block of particle scatter counts to the specified output file.

This function writes a chunk of integer data, corresponding to block_size timesteps, to the given binary file. Each integer represents the number of times a particle scattered in a given timestep.

Parameters:

• filename -- [in] Path to the output binary scatter count data file.

• npts -- [in] Number of particles.

• block_size -- [in] Number of timesteps of data contained in scat_count_block.

• scat_count_block -- [in] Pointer to the block of scatter count data (int array).

static void calculate_system_energies(double **particles, int npts, double deltaM, double *total_KE_out, double *total_PE_out)

Calculates the total kinetic and potential energy of the N-body system.

This function iterates over all particles to compute the system's total kinetic energy (KE) and potential energy (PE) at a specific snapshot in time. The particles array must be sorted by radius before calling this function to ensure the rank-based potential energy calculation is correct.

Parameters:

• particles -- [in] The main particle data array, sorted by radius.

• npts -- [in] The total number of particles in the simulation.

• deltaM -- [in] The mass of a single particle ( \(M_{\odot}\)).

• total_KE_out -- [out] Pointer to store the calculated total kinetic energy.

• total_PE_out -- [out] Pointer to store the calculated total potential energy.

static int check_and_warn_negative_fQ(gsl_interp **fofE_interp_ptr, gsl_interp_accel **fofE_acc_ptr, double *E_array, double *I_array, int *n_points_ptr, const char *label, const char *profile_name, int verbose)

Checks for negative distribution function values and corrects minor artifacts.

Evaluates \(f(E) = dI/dE\) at each point in the original arrays. If negative \(f(E)\) values are found (where \(I(E)\) decreases):

• For < 20 negative points: Removes them, rebuilds spline, warns (verbose mode only)

• For \(\geq 20\) negative points: Prompts user to abort or continue (verbose mode only)

Parameters:

• fofE_interp_ptr -- Pointer to GSL interpolator (rebuilt if correction applied)

• fofE_acc_ptr -- Pointer to GSL accelerator (rebuilt if correction applied)

• E_array -- Array of energy values (modified in-place during correction)

• I_array -- Array of \(I(E)\) values (modified in-place during correction)

• n_points_ptr -- Pointer to number of points (updated during correction)

• label -- String label "f(E)" or "f(Q)" for display

• profile_name -- String name of profile (e.g., "NFW", "Cored")

• verbose -- If 1: show warnings/prompts. If 0: silent correction (diagnostic mode)

Returns:

int 1 to abort, 0 to continue

static int check_strict_monotonicity(const double *arr, int n, const char *name)

Check if an array is strictly monotonically increasing.

Parameters:

• arr -- Array to check

• n -- Number of elements

• name -- Name of the array for debug messages

Returns:

1 if strictly monotonic, 0 otherwise

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.

static void cleanup_trajectory_buffers(int num_traj_particles, int nlowest)

Frees all memory associated with the trajectory buffering system.

This function is called at the end of the simulation. It ensures any remaining data in the buffers is flushed, closes all open file handles, and deallocates all buffer memory.

Parameters:

• num_traj_particles -- Number of low-ID particles tracked.

• nlowest -- Number of low-L particles tracked.

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

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

Sorts LAndIndex structures in ascending order based on their 'L' (angular momentum or squared difference from a reference \(L\)) value. Used with qsort for ordering particles by their \(L\) values, typically for selecting particles with lowest \(L\) or \(L\) closest to a target.

Parameters:

• a -- [in] Pointer to the first LAndIndex structure.

• b -- [in] Pointer to the second LAndIndex structure.

Returns:

int -1 if a->L < b->L, 1 if a->L > b->L, 0 if equal.

int compare_pair_by_first_element(const void *a, const void *b)

Comparison function for qsort to sort RrPsiPair structures by the 'rr' member.

Used to sort an array of RrPsiPair structures in ascending order based on their radial (rr) component. This is primarily used when preparing data for splines where the x-axis (e.g., radius or -Psi) must be strictly monotonic.

Parameters:

• a -- Pointer to the first RrPsiPair structure.

• b -- Pointer to the second RrPsiPair structure.

Returns:

int -1 if pa->rr < pb->rr, 1 if pa->rr > pb->rr, 0 otherwise.

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

Comparison function for qsort to sort PartData structures by radius.

Compares two PartData structures based on their rad (radius) member for sorting in ascending order. Handles NaN values by placing them consistently at the beginning (NaN values sorted first).

Parameters:

• a -- Pointer to the first PartData structure.

• b -- Pointer to the second PartData structure.

Returns:

int -1 if pa->rad < pb->rad, 1 if pa->rad > pb->rad, 0 otherwise. Specific return for NaNs ensures consistent ordering.

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

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

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

Particle data structure: columns[particle_index][component_index]. Comparison based on columns[i][0] (radial position).

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).

double cored_potential_wrapper(double r, void *params)

Wrapper function for Cored Plummer-like potential calculation.

Parameters:

• r -- Radius (kpc).

• params -- Void pointer to cored_potential_params struct.

Returns:

Potential value at radius r.

static int correct_negative_fE_and_rebuild_IE(double *fE_values, int *is_negative, double *E_array, double *I_array, int n_points, const char *label)

Corrects negative \(f(E)\) regions by linear interpolation and reconstructs \(I(E)\).

Given \(f(E)\) values with some negative regions, linearly interpolates \(f(E)\) over consecutive negative ranges, then integrates to reconstruct \(I(E)\).

Parameters:

• fE_values -- Array of \(f(E)\) values (modified in-place)

• is_negative -- Array marking negative \(f(E)\) points (1 = negative, 0 = positive)

• E_array -- Array of energy values

• I_array -- Array to store reconstructed \(I(E)\) values (output)

• n_points -- Number of points

• label -- Label for logging ("f(E)" or "f(Q)")

Returns:

Number of ranges corrected

static int create_backup_file(const char *source_file)

Creates a backup of a file with a .backup extension.

Creates a byte-for-byte copy of the source file, preserving the all_particle_data file state before truncation during restart operations.

Parameters:

source_file -- Path to the source file to back up.

Returns:

Returns 0 on success, -1 on failure.

threevector crossproduct(threevector X, threevector Y)

Computes the vector cross product of two three-dimensional vectors.

Calculates \(Z = X \times Y\), where \(X = (X_x, X_y, X_z)\) and \(Y = (Y_x, Y_y, Y_z)\). The components of the resulting vector \(Z\) are determined by: \(Z_x = X_y Y_z - X_z Y_y\) \(Z_y = X_z Y_x - X_x Y_z\) \(Z_z = X_x Y_y - X_y Y_x\) This follows the standard right-hand rule for vector cross products.

Parameters:

• X -- [in] The first threevector operand.

• Y -- [in] The second threevector operand.

Returns:

threevector The resulting vector \(Z = X \times Y\).

double density_derivative_om_cored(double r)

Osipkov-Merritt augmented density functions for anisotropic models.

The OM model uses radially-varying anisotropy \(\beta(r) = r^2/(r^2 + r_a^2)\) implemented via augmented density \(\rho_Q(r) = \rho(r)(1 + r^2/r_a^2)\). These functions compute derivatives of the augmented density for use in Eddington-like inversion to obtain the distribution function \(f(Q)\).

Calculates \(d\rho_Q/dr\) for the Osipkov-Merritt Cored Plummer-like profile.

Computes the radial derivative of the OM augmented density \(\rho_Q(r) = \rho(r)(1 + r^2/r_a^2)\). Uses the product rule: \(d\rho_Q/dr = (d\rho/dr)(1+r^2/r_a^2) + \rho(r)(2r/r_a^2)\). This is used in the Eddington-like inversion to find \(f(Q)\).

Parameters:

r -- [in] Radial coordinate (kpc).

Returns:

double The value of \(d\rho_Q/dr\) at radius r.

double density_derivative_om_nfw(double r, double rc_param, double nt_nfw_scaler, double falloff_C_param)

Calculates \(d\rho_Q/dr\) for the Osipkov-Merritt NFW-like profile.

Computes the radial derivative of the OM augmented density \(\rho_Q(r) = \rho(r)(1 + r^2/r_a^2)\). Uses the product rule: \(d\rho_Q/dr = (d\rho/dr)(1+r^2/r_a^2) + \rho(r)(2r/r_a^2)\). This is used in the Eddington-like inversion to find \(f(Q)\).

Parameters:

• r -- [in] Radial coordinate (kpc).

• rc_param -- [in] Scale radius \(r_c\) of the NFW-like profile (kpc).

• nt_nfw_scaler -- [in] Density normalization constant (nt_nfw) for the profile.

• falloff_C_param -- [in] Falloff transition factor C for the power-law cutoff.

Returns:

double The value of \(d\rho_Q/dr\) at radius r.

static inline double density_hernquist(double r, double M, double a)

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

Analytical density for the Hernquist profile.

Orchestrates the overall simulation workflow:

1. Parses command-line arguments.

2. Sets up global parameters and logging.

3. Initializes random number generators.

4. Generates or loads initial conditions (ICs) for NFW, Cored Plummer-like, or Hernquist profiles.

   • Includes theoretical calculations for density, mass, potential, and \(f(E)\) splines.

   • Includes a diagnostic loop to test IC generation with varied numerical parameters.

   • Performs particle sampling based on the derived distribution function.

5. Optionally performs tidal stripping and re-assigns particle IDs.

6. Converts particle velocities to physical simulation units.

7. Executes the main N-body timestepping loop using a selected integration method.

   • Performs gravitational updates.

   • Optionally performs SIDM scattering via handle_sidm_step.

   • Tracks particle trajectories and energies.

   • Periodically writes simulation data and progress.

8. If g_doAllParticleData is enabled, processes all particle data to generate snapshot files for Rank/Mass/Radius/Velocity/Potential/Energy/Density.

9. Writes final summary plots and theoretical profiles.

10. Cleans up allocated resources. Handles restart/resume functionality by checking for existing data products.

Note

This application supports OpenMP for parallelization in various sections.

Warning

Large particle counts or long simulations can be memory and CPU intensive. Disk space requirements for full data output can also be significant.

Parameters:

• argc -- [in] Standard argument count from the command line.

• argv -- [in] Standard array of argument strings from the command line.

• r -- Radius (kpc).

• M -- Total halo mass ( \(M_{\odot}\)).

• a -- Scale radius \(a\) of the Hernquist profile (kpc).

Returns:

int Exit code: 0 for successful execution, non-zero for errors.

Returns:

Physical density ( \(M_{\odot}\) / kpc \(^3\)).

static inline double df_hernquist_aniso(double E_bind, double L, double M, double a)

Anisotropic distribution function for Hernquist profile with arbitrary beta.

Parameters:

• E_bind -- Binding energy, E_bind = -E_total (must be positive). Units are (kpc/Myr) \(^2\).

• L -- Angular momentum magnitude (kpc \(^2\)/Myr).

• M -- Total halo mass ( \(M_{\odot}\)).

• a -- Scale radius \(a\) of the Hernquist profile (kpc).

Returns:

Value of the distribution function f(E, L).

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/(\sigma\sqrt{2\pi})) \exp(-0.5(x/\sigma)^2)\), with x being the distance log_r_grid[i] - log_r_grid[j] and \(\sigma\) = 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^2)\)) which makes it slower for large grid_size.

See also

fft_gaussian_convolution

See also

gaussian_convolution

Note

Resulting smoothed density is clamped to a minimum value of 1e-10. This function is inherently thread-safe as it only reads inputs and writes to distinct elements of the output array without shared intermediate state.

Parameters:

• density_grid -- [in] Input density grid array (values corresponding to log_r_grid).

• grid_size -- [in] Number of points in the input grid and density arrays.

• log_r_grid -- [in] Array of logarithmic radial grid coordinates (log10(r)). Can be non-uniformly spaced.

• sigma_log -- [in] Width (standard deviation) of the Gaussian kernel in log10-space.

• result -- [out] Output array (pre-allocated, size grid_size) for the smoothed density field.

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 over a physical timestep h.

This function implements an adaptive leapfrog algorithm to advance a particle's state \((r, v_{rad})\) over a full physical timestep h ( \(\Delta T_{phys}\)). It iteratively refines the integration by comparing a "coarse" integration (using \(2N+1\) micro-steps via doMicroLeapfrog) with a "fine" integration (using \(4N+1\) micro-steps). The subdivision factor \(N\) starts at 1 and is doubled if the relative differences in final radius and velocity between coarse and fine passes exceed radius_tol and velocity_tol, respectively. This process repeats up to a maximum subdivision factor max_subdiv. The final state \((r_{out}, v_{out})\) for the step h is chosen based on out_type (coarse, fine, or Richardson extrapolation) once convergence is met or max_subdiv is reached.

Parameters:

• i -- [in] Particle index (0 to npts-1), for rank in gravitational force calculation.

• npts -- [in] Total number of particles in the simulation.

• r_in -- [in] Initial radial position (kpc) at the start of the step h.

• v_in -- [in] Initial radial velocity (kpc/Myr) at the start of h.

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr), conserved.

• h -- [in] Full physical timestep size \(\Delta T_{phys}\) (Myr).

• radius_tol -- [in] Relative convergence tolerance for radius comparison.

• velocity_tol -- [in] Relative convergence tolerance for velocity comparison.

• max_subdiv -- [in] Maximum allowed subdivision factor \(N\) for micro-steps.

• grav -- [in] Gravitational constant G (simulation units, e.g., G_CONST).

• out_type -- [in] Result selection mode for converged integration: 0 for coarse result, 1 for fine result, 2 for Richardson extrapolation.

• r_out -- [out] Pointer to store the final radial position (kpc) after time h.

• v_out -- [out] Pointer to store the final radial velocity (kpc/Myr) after time h.

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 a full adaptive leapfrog step using Levi-Civita regularization over a physical time interval dt.

This function integrates a particle's motion over a physical timestep dt ( \(\Delta T_{phys}\)) by taking multiple adaptive steps in fictitious Levi-Civita time \(τ\). It repeatedly calls doSingleTauStepAdaptiveLeviCivita to advance the state in \((ρ, v_{rad}, t_{phys})\) coordinates, where \(ρ = \sqrt{r}\). Each call to doSingleTauStepAdaptiveLeviCivita takes an adaptive \(Δτ\) step. The loop continues until the accumulated physical time t_cur (from summing \(Δt_{phys}\) corresponding to each \(Δτ\)) reaches or exceeds the target dt. If a \(Δτ\) step overshoots dt, linear interpolation is used to find the state precisely at t_cur = dt. The final regularized state \((ρ_f, v_f)\) is then transformed back to physical coordinates \((r_{out}, v_{out})\). An initial guess for \(Δτ\) is made based on N_taumin and the initial radius.

Parameters:

• i -- [in] Particle index (0 to npts-1), used for rank in gravitational force calculation.

• npts -- [in] Total number of particles in the simulation.

• r_in -- [in] Initial physical radial position (kpc) at the start of the physical step dt.

• v_in -- [in] Initial physical radial velocity (kpc/Myr) at the start of dt.

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr), conserved during integration.

• dt -- [in] The full physical timestep \(\Delta T_{phys}\) (Myr) to advance the particle.

• N_taumin -- [in] Target number of fictitious \(τ\)-steps within dt; influences the initial \(Δτ\) guess.

• radius_tol -- [in] Relative convergence tolerance for \(ρ\) comparison within each adaptive \(τ\)-step.

• velocity_tol -- [in] Relative convergence tolerance for velocity comparison within each adaptive \(τ\)-step.

• max_subdiv -- [in] Maximum allowed subdivision factor N for micro-steps within each adaptive \(τ\)-step.

• grav -- [in] Gravitational constant G (simulation units, e.g., G_CONST).

• out_type -- [in] Result selection mode for micro-steps within doSingleTauStepAdaptiveLeviCivita: 0 for coarse, 1 for fine, 2 for Richardson extrapolation.

• r_out -- [out] Pointer to store the final physical radial position (kpc) after time dt.

• v_out -- [out] Pointer to store the final physical radial velocity (kpc/Myr) after time dt.

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 over a physical time dt using Levi-Civita regularization.

This function implements a leapfrog-like integration scheme in Levi-Civita regularized coordinates \((ρ, v_{rad})\) and fictitious time \(τ\). The physical coordinates are \(r = ρ^2\), and physical time \(t_{phys}\) is related to \(τ\) by \(dt_{phys} = ρ^2 dτ\). The integration proceeds by taking variable \(Δτ\) steps (estimated based on N_taumin) until the accumulated physical time t_phys reaches or exceeds the target physical timestep dt. If a step overshoots dt, linear interpolation is used to obtain the state precisely at the target physical time. This method is particularly effective for handling close encounters where \(r → 0\).

Parameters:

• i -- [in] Particle index (0 to npts-1), for rank in force calculation.

• npts -- [in] Total number of particles.

• r_in -- [in] Initial physical radial position (kpc) at the start of the physical step dt.

• v_in -- [in] Initial physical radial velocity (kpc/Myr) at the start of dt.

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr) for the force calculation.

• dt -- [in] The full physical timestep \(\Delta T_{phys}\) (Myr) to advance the particle.

• N_taumin -- [in] Target number of fictitious \(τ\)-steps within the dt interval; influences the initial guess for \(Δτ\).

• grav -- [in] Gravitational constant G (simulation units).

• r_out -- [out] Pointer to store the final physical radial position (kpc) after time dt.

• v_out -- [out] Pointer to store the final physical radial velocity (kpc/Myr) after time dt.

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 leapfrog integration step using a fixed number of micro-steps.

This helper function implements the leapfrog (Kick-Drift-Kick) integration over a total time interval h by dividing it into a sequence of micro-steps. The number of micro-steps is subSteps (e.g., \(2N+1\) for a "coarse" pass or \(4N+1\) for a "fine" pass in an adaptive scheme, where \(N\) is N_subdivision_factor). The micro-timestep sizes for kicks and drifts are adjusted based on N_subdivision_factor and whether it is a coarse or fine integration sequence. Sequence: Initial half-kick, \(((\text{subSteps}-1)/2 - 1)\) full Drift-Kick pairs, a final full Drift, and a final half-Kick.

Parameters:

• i -- [in] Particle index (0 to npts-1), for rank in gravitational force.

• npts -- [in] Total number of particles.

• r_in -- [in] Input radial position (kpc) at the start of the total interval h.

• v_in -- [in] Input radial velocity (kpc/Myr) at the start of h.

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr).

• h -- [in] Total physical time interval for this leapfrog sequence (Myr).

• N_subdivision_factor -- [in] Base subdivision factor \(N\) used to determine micro-timestep sizes.

• subSteps -- [in] Total number of Kicks/Drifts (e.g., \(2N+1\) or \(4N+1\)).

• grav -- [in] Gravitational constant G (simulation units).

• r_out -- [out] Pointer to store the output radial position (kpc) after time h.

• v_out -- [out] Pointer to store the output radial velocity (kpc/Myr) after time h.

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 using a fixed number of micro-steps.

This helper function advances the particle state in regularized coordinates \((ρ, v_{rad}, t_{phys})\) over a total fictitious time interval h_tau ( \(Δτ_{total}\)) by taking a specified number of subSteps fixed-size micro-steps ( \(δτ = Δτ_{total} / \text{subSteps}\)). Each micro-step uses a leapfrog-like scheme (Kick-Drift-Kick for \(ρ, v_{rad}\)) and updates the accumulated physical time \(t_{phys}\) using \(dt_{phys} = δτ \cdot ρ_{mid}^2\). This function is called by doSingleTauStepAdaptiveLeviCivita for its coarse and fine passes.

Parameters:

• i -- [in] Particle index (0 to npts-1), for rank in force calculation.

• npts -- [in] Total number of particles.

• rho_in -- [in] Initial \(ρ = \sqrt{r}\) at the start of the h_tau interval.

• v_in -- [in] Initial radial velocity \(v_{rad}\) at the start of h_tau.

• t_in -- [in] Initial accumulated physical time \(t_{phys}\) at the start of h_tau.

• subSteps -- [in] Number of fixed micro-steps to perform over h_tau.

• h_tau -- [in] Total fictitious time interval \(Δτ_{total}\) for this integration sequence.

• grav -- [in] Gravitational constant G (simulation units).

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr).

• rho_out -- [out] Pointer to store the final \(ρ\) after h_tau.

• v_out -- [out] Pointer to store the final \(v_{rad}\) after h_tau.

• t_out -- [out] Pointer to store the final accumulated \(t_{phys}\) after h_tau.

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 over a proposed fictitious time h_guess.

This function integrates the equations of motion in regularized \((ρ, v_{rad}, t_{phys})\) coordinates over a proposed fictitious time interval h_guess ( \(Δτ_{guess}\)). It employs an adaptive refinement strategy by comparing a "coarse" integration (using \(2N+1\) micro-steps via doMicroLeviCivita) with a "fine" integration (using \(4N+1\) micro-steps). The subdivision factor \(N\) starts at 1 and is doubled if the relative differences in \(ρ\) and \(v_{rad}\) between coarse and fine results exceed radius_tol and velocity_tol, respectively. This continues up to max_subdiv. The final state for the \(Δτ_{guess}\) step is chosen based on out_type (coarse, fine, or Richardson extrapolation). The function outputs the final \(ρ_{out}\), \(v_{out}\) (radial velocity), and accumulated physical time \(t_{out}\) corresponding to this adaptive \(Δτ\) step.

Parameters:

• i -- [in] Particle index (0 to npts-1), for rank in force calculation.

• npts -- [in] Total number of particles.

• rho_in -- [in] Initial regularized radial coordinate \(ρ = \sqrt{r}\) at the start of \(Δτ_{guess}\).

• v_in -- [in] Initial radial velocity \(v_{rad}\) at the start of \(Δτ_{guess}\).

• t_in -- [in] Initial accumulated physical time \(t_{phys}\) at the start of \(Δτ_{guess}\).

• h_guess -- [in] The proposed total fictitious time step \(Δτ_{guess}\) for this adaptive step.

• radius_tol -- [in] Relative convergence tolerance for \(ρ\) comparison.

• velocity_tol -- [in] Relative convergence tolerance for \(v_{rad}\) comparison.

• max_subdiv -- [in] Maximum subdivision factor \(N\) for micro-steps within doMicroLeviCivita.

• grav -- [in] Gravitational constant G (simulation units).

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr).

• out_type -- [in] Result selection for converged micro-integration: 0=coarse, 1=fine, 2=Richardson.

• rho_out -- [out] Pointer to store the final \(ρ\) after the adaptive \(Δτ_{guess}\) step.

• v_out -- [out] Pointer to store the final \(v_{rad}\) after the adaptive \(Δτ_{guess}\) step.

• t_out -- [out] Pointer to store the final accumulated physical time \(t_{phys}\) after this step.

double dotproduct(threevector X, threevector Y)

Computes the scalar dot product of two three-dimensional vectors.

Calculates \(X \cdot Y = X_x Y_x + X_y Y_y + X_z Y_z\). The dot product is a measure of the projection of one vector onto another and is used in various physics calculations, such as determining the magnitude squared of a vector ( \(V \cdot V = |V|^2\)) or the angle between vectors.

Parameters:

• X -- [in] The first threevector operand.

• Y -- [in] The second threevector operand.

Returns:

double The scalar result of the dot product \(X \cdot Y\).

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

Comparison function for sorting double values in ascending order.

This function is designed to be used with qsort or other standard library sorting functions that require a comparator. It takes two void pointers, casts them to const double*, dereferences them, and compares their values.

Parameters:

• a -- [in] Pointer to the first double value.

• b -- [in] Pointer to the second double value.

Returns:

int -1 if the first double is less than the second, 1 if the first double is greater than the second, 0 if they are equal.

double drhodr(double r)

Calculates \(d\rho/dr\) for the Cored Plummer-like density profile.

The density profile is \(\rho(r) \propto (1 + (r/r_c)^2)^{-3}\). This function computes its analytical derivative with respect to \(r\). It directly uses the RC macro for the scale radius.

Parameters:

r -- [in] Radial coordinate (kpc) at which to evaluate the derivative.

Returns:

double The value of \(d\rho/dr\) at radius r.

double drhodr_hernquist(double r, double a_scale, double M_total)

Calculates \(d\rho/dr\) for the Hernquist profile.

Computes the radial derivative of the Hernquist density profile \(\rho(r) = M_{\text{total}} a / [2\pi r (r+a)^3]\). The derivative is: \(d\rho/dr = -(M a / 2\pi) (4r + a) / [r^2(r+a)^4]\)

Parameters:

• r -- [in] Radial coordinate (kpc).

• a_scale -- [in] Scale radius \(a\) of the Hernquist profile (kpc).

• M_total -- [in] Total mass normalization (not used in shape calculation).

Returns:

double The value of \(d\rho/dr\) at radius r.

double drhodr_om_hernquist(double r, double a_scale, double M_total)

Calculates \(d\rho_Q/dr\) for the Osipkov-Merritt Hernquist profile.

Computes the radial derivative of the OM augmented density \(\rho_Q(r) = \rho(r)(1 + r^2/r_a^2)\) for the Hernquist profile. Uses the product rule: \(d\rho_Q/dr = (d\rho/dr)(1+r^2/r_a^2) + \rho(r)(2r/r_a^2)\).

Parameters:

• r -- [in] Radial coordinate (kpc).

• a_scale -- [in] Scale radius \(a\) of the Hernquist profile (kpc).

• M_total -- [in] Total mass normalization (not used in shape calculation).

Returns:

double The value of \(d\rho_Q/dr\) at radius r.

double drhodr_profile_nfwcutoff(double r, double rc_param, double nt_nfw_scaler, double falloff_C_param)

Calculates \(d\rho/dr\) for the NFW-like profile with a power-law cutoff.

The NFW-like density profile used is: \(\rho(r) = \text{nt_nfw_scaler} \times [ (r_s + \epsilon)(1+r_s)^2 (1 + (r_s/C)^N) ]^{-1}\) where \(r_s = r / \text{rc_param}\), \(\epsilon\) is a softening parameter (0.01), \(C\) is the falloff_C_param, and \(N\) is a power-law index (10.0). This function computes the analytical derivative of this \(\rho(r)\) with respect to \(r\).

Parameters:

• r -- [in] Radial coordinate (kpc) at which to evaluate the derivative.

• rc_param -- [in] Scale radius \(r_c\) of the NFW-like profile (kpc).

• nt_nfw_scaler -- [in] Density normalization constant (nt_nfw) for the profile.

• falloff_C_param -- [in] Falloff transition factor \(C\) for the power-law cutoff.

Returns:

double The value of \(d\rho/dr\) at radius r. Returns 0.0 if rc_param is non-positive.

static inline double dRhoDtaufun(double rhoVal, double vVal)

Calculates \(d\rho/d\tau\), the derivative of the regularized coordinate \(\rho\) with respect to fictitious time \(\tau\).

In Levi-Civita regularization, \(d\rho/d\tau = (1/2) \rho v_{rad}\), where \(\rho = \sqrt{r}\) and \(v_{rad}\) is the radial velocity in physical units (though often represented as \(v\) or \(v_{\rho}\) in transformed equations of motion depending on the specific formulation). This function implements this relationship.

Parameters:

• rhoVal -- [in] The current value of the regularized radial coordinate \(\rho = \sqrt{r}\).

• vVal -- [in] The current radial velocity \(v_{rad}\) (kpc/Myr).

Returns:

double The value of \(d\rho/d\tau\) (kpc \(^{3/2}\)/Myr).

static inline double effective_angular_force(double r, double ell)

Computes the effective centrifugal acceleration due to angular momentum.

Parameters:

• r -- [in] Radius (kpc).

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr).

Returns:

Centrifugal acceleration: \(L^2/r^3\) (kpc/Myr \(^2\)).

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.

See also

effective_angular_force

See also

doLeviCivitaLeapfrog

Parameters:

• rho -- [in] Transformed radial coordinate (sqrt(r)). Units: sqrt(kpc).

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr).

Returns:

\(dv/d\tau\) centrifugal term in Levi-Civita coordinates (kpc \(^2\)/Myr \(^2\)).

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.

Note

Calls cleanup_all_particle_data() before exiting to free allocated memory.

Warning

This function does not return; execution is terminated.

Parameters:

• msg -- [in] The error message to display.

• arg -- [in] Optional argument value that caused the error (shown in quotes), or NULL to omit.

• prog -- [in] The program name to display in the &#8212;help usage suggestion.

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

Safely evaluates a GSL spline with robust bounds checking.

Evaluates the spline at the given value. If the value is outside the spline's defined domain, this function clamps it to the nearest valid boundary to prevent a GSL domain error. Handles NULL pointers for the spline or accelerator by logging an error and returning 0.0.

Parameters:

• spline -- Pointer to the initialized GSL spline object.

• acc -- Pointer to the GSL interpolation accelerator.

• value -- The value at which to evaluate the spline.

Returns:

The interpolated value from the spline. If value is out of bounds, returns the value at the nearest boundary. If spline or acc is NULL, prints an error to stderr and returns 0.0.

static inline double evaluatespline_with_boundary(gsl_spline *spline, gsl_interp_accel *acc, double x, double x_min, double x_max, double y_min, double y_max)

Evaluates a spline, returning boundary values when outside domain.

Parameters:

• spline -- GSL spline object

• acc -- GSL interpolation accelerator

• x -- Point to evaluate

• x_min -- Minimum x value in spline

• x_max -- Maximum x value in spline

• y_min -- Value to return when x < x_min

• y_max -- Value to return when x > x_max

Returns:

Interpolated value or boundary value

double fEintegrand(double t, void *params)

Calculates the integrand for the distribution function calculation.

Uses splines to evaluate the potential and mass at a given energy and radius, and computes the distribution function integrand according to Eddington's formula.

Parameters:

• t -- Integration variable

• params -- Structure containing necessary splines and energy value

Returns:

Value of the distribution function integrand at the specified point

double fEintegrand_hernquist(double t_integration_var, void *params)

Integrand for the NFW distribution function \(I(E)\) calculation.

The function computes the integrand \(2 \cdot (d\rho/d\Psi)\) for Eddington's formula, using the transformed integration variable \(t = \sqrt{E_{shell} - \Psi(r)}\). It determines the radius r corresponding to a given t by first finding \(\Psi(r) = E_{shell} - t^2\) and then using a spline for \(r(\Psi)\). The derivatives \(d\rho/dr\) and \(d\Psi/dr\) are then calculated at that radius.

Integrand of Hernquist-specific Eddington inversion for computing \(I(E)\) or \(I(Q)\).

The integrand is \(2(d\rho/d\Psi)\) which is positive due to sign convention \(d\Psi/dr = -GM/r^2\). Uses transformed variable t_integration_var = sqrt(E_shell - Psi_true) to avoid singularities. For OM model, E_shell represents \(Q = E - L^2/(2r_a^2)\).

Parameters:

• t_integration_var -- The integration variable, \(t = \sqrt{E_{shell} - \Psi(r)}\).

• params -- Pointer to a fE_integrand_params_NFW_t struct containing the energy shell, splines, and profile parameters.

• t_integration_var -- [in] The t variable where t = sqrt(E_shell - Psi_true). Integration bounds are [0, sqrt(E_shell - Psimin_global)].

• params -- [in] Pointer to fE_integrand_params_hernquist_t containing E_shell, splines, profile parameters, and physical range limits.

Returns:

The value of the integrand \(2 \cdot (d\rho/d\Psi)\). Returns 0.0 if the potential is outside the physical range, the radius is non-physical, \(d\Psi/dr\) is near zero, or the result is non-finite.

Returns:

double The value of the integrand \(2(d\rho/d\Psi)\). Returns 0.0 if invalid.

double fEintegrand_nfw(double t_integration_var, void *params)

Integrand for NFW distribution function \(I(E)\) calculation.

Parameters:

• t_integration_var -- Integration variable \(t = \sqrt{E_{shell} - \Psi(r)}\).

• params -- Pointer to fE_integrand_params_NFW_t struct.

Returns:

Integrand value \(2(d\rho/d\Psi)\).

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).

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.

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.

Parameters:

• density_grid -- [in] Input density grid array (values corresponding to log_r_grid).

• grid_size -- [in] Number of points in the input grid and density arrays.

• log_r_grid -- [in] Array of logarithmic radial grid coordinates (log10(r)). Must be uniformly spaced.

• sigma_log -- [in] Width (standard deviation) of the Gaussian kernel in log10-space.

• result -- [out] Output array (pre-allocated, size grid_size) for the smoothed density field.

static void finalize_debug_energy_output(void)

Writes collected debug energy data to debug_energy_compare.dat.

Outputs time evolution of theoretical vs dynamical energy for the tracked particle (DEBUG_PARTICLE_ID), including KE and PE components. Called at simulation end if g_doDebug is enabled.

static void finalize_energy_diagnostics(int noutsnaps)

Writes the collected total system energy diagnostics to a file.

This function is called at the end of the simulation. It iterates through the global snapshot arrays and writes the time evolution of the system's total kinetic, potential, and combined energy to a file named data/total_energy_vs_time<suffix>.dat.

Note

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

Parameters:

noutsnaps -- [in] The total number of snapshots recorded.

static int find_last_common_complete_snapshot(int npts, int dtwrite, int num_traj_particles, int nlowest, const char *file_suffix)

Find the last complete snapshot that exists across ALL output files.

When a simulation is interrupted, different files may have different amounts of data due to buffering. This function finds the last snapshot that is complete in ALL files to ensure consistency on restart.

Parameters:

• npts -- Number of particles

• dtwrite -- Timestep interval for writes

• num_traj_particles -- Number of trajectory particles

• nlowest -- Number of lowest-L particles

• file_suffix -- The file suffix string

Returns:

The last complete snapshot number common to all files, or -1 on error

static int find_last_processed_snapshot(int *snapshot_steps, int noutsnaps)

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

This function is called when --restart is active. It checks for the existence and basic integrity of snapshot data files (e.g., Rank_Mass_Rad_VRad_unsorted_t%05d.dat and Rank_Mass_Rad_VRad_sorted_t%05d.dat) corresponding to the timesteps listed in snapshot_steps. Integrity is checked by comparing file sizes against those of the first snapshot's files (snapshot_steps[0]), allowing for a small percentage tolerance (typically +/- 5%). The function aims to determine from which snapshot index s (in snapshot_steps) the post-processing (e.g., Rank file generation) should resume. It assumes g_file_suffix is correctly set to identify the relevant run's files.

Parameters:

• snapshot_steps -- [in] Array of integer timestep numbers for which snapshot data was expected to be written (often these are the dtwrite intervals).

• noutsnaps -- [in] The total number of snapshot indices in the snapshot_steps array.

Returns:

int The index s into snapshot_steps corresponding to the first snapshot that needs to be processed (i.e., last valid snapshot index + 1).

• Returns -1 if no valid snapshot files are found (implying processing should start from index 0).

• Returns -2 if all expected snapshot files for all unique snapshot numbers exist and appear valid (implying no further snapshot processing is needed for this phase). Prints status messages to stdout and logs warnings/errors.

double find_minimum_useful_radius(double (*potential_func)(double, void*), void *params, double scale_radius, double tolerance)

Finds the minimum useful radius where potential values are numerically distinguishable.

Starting from a safe radius, steps down logarithmically towards \(r=0\) until consecutive potential values differ by less than the tolerance.

Parameters:

• potential_func -- Function pointer to evaluate potential at radius \(r\)

• params -- Parameters to pass to potential_func

• scale_radius -- Characteristic scale of the profile

• tolerance -- Relative tolerance for considering potentials equal

Returns:

Minimum radius where potential is numerically distinguishable from smaller radii

static void flush_double_buffer_to_disk(const char *tag)

Flushes double-precision buffer to disk files.

Writes the entire 4-snapshot buffer to disk, overwriting previous buffer files. Creates two binary files:

• double_buffer_all_particle_data_<tag>.dat (Rank, R, Vrad, L)

• double_buffer_all_particle_phi_<tag>.dat (phi)

File format per particle per snapshot:

• data file: int32 Rank + double R + double Vrad + double L = 28 bytes

• phi file: double phi = 8 bytes

Parameters:

tag -- Simulation tag for filename

static void flush_trajectory_buffers(int items_to_write, int num_traj_particles, int nlowest)

Flushes all buffered trajectory data to their respective files.

This function writes the collected trajectory data for a block of timesteps to disk. It handles opening files on the first write and appends data on subsequent calls. It is synchronized with the main snapshot buffer flush.

Parameters:

• items_to_write -- The number of timesteps currently in the buffer to write.

• num_traj_particles -- The number of low-ID particles being tracked.

• nlowest -- The number of low-L particles being tracked.

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

Calculates the total effective force per unit mass in Levi-Civita transformed coordinates.

This function computes \(F_{\rho}/m = (F_{grav,\rho} + F_{centrifugal,\rho})/m\), where \(F_{grav,\rho}\) is the gravitational force and \(F_{centrifugal,\rho}\) is the effective centrifugal force, both expressed in the regularized radial coordinate \(\rho = \sqrt{r}\). It calls gravitational_force_rho_v and effective_angular_force_rho_v. This combined force is used in the equations of motion for Levi-Civita regularization.

Parameters:

• i -- [in] Particle index (0 to npts-1), for rank in gravitational force calculation.

• npts -- [in] Total number of particles.

• totalmass -- [in] Total halo mass of the system ( \(M_{\odot}\)) used for gravitational force.

• grav -- [in] Gravitational constant G (simulation units).

• ell -- [in] Angular momentum per unit mass (kpc \(^2\)/Myr).

• rhoVal -- [in] Current value of the regularized radial coordinate \(\rho = \sqrt{r}\).

Returns:

double \(dv/d\tau\) in Levi-Civita coordinates (kpc \(^2\)/Myr \(^2\)).

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.

Parameters:

• size_in_bytes -- [in] The size in bytes to format.

• buffer -- [out] Output buffer for the formatted string.

• buffer_size -- [in] Size of the output buffer.

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.

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.

See also

fscanf_bin

Note

Skips optional width/precision specifiers in the format string.

Parameters:

• fp -- [in] File pointer to write to (must be opened in binary mode).

• format -- [in] Format string with specifiers (d, f, g, e). Other characters are ignored.

• ... -- [in] Variable arguments matching the format specifiers.

Returns:

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

int fprintf_bin_dbl(FILE *fp, const char *format, ...)

Writes binary data to a file using a printf-like format string (DOUBLE PRECISION).

Writes floating-point values as 8-byte doubles, preserving full precision. Avoids the precision loss from double→float32→double conversion.

Used exclusively for the double_buffer_all_particle_* files which maintain the last 4 snapshots in full precision for accurate restart/extend operations.

See also

fprintf_bin (float32 version)

See also

fscanf_bin_dbl

Note

Integer types (d) written as int (4 bytes), floats (f/g/e) written as double (8 bytes).

Parameters:

• fp -- [in] File pointer to write to (must be opened in binary mode).

• format -- [in] Format string with specifiers (d, f, g, e). Other characters are ignored.

• ... -- [in] Variable arguments matching the format specifiers.

Returns:

The number of items written according to 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.

See also

fprintf_bin

Note

Skips optional width/precision specifiers in the format string.

Parameters:

• fp -- [in] File pointer to read from (must be opened in binary mode).

• format -- [in] Format string with specifiers (d, f, g, e). Other characters are ignored.

• ... -- [out] Variable pointer arguments matching the format specifiers (e.g., int*, double*).

Returns:

The number of items successfully read and assigned. Stops reading on first failure or EOF.

int fscanf_bin_dbl(FILE *fp, const char *format, ...)

Reads binary data from a file using a scanf-like format string (DOUBLE PRECISION).

Reads floating-point values as 8-byte doubles, preserving full precision. Avoids the precision loss from float32→double conversion when loading snapshot data.

Used to restore full-precision snapshot data during restart/extend operations, avoiding the precision loss inherent in float32 storage.

See also

fscanf_bin (float32 version)

See also

fprintf_bin_dbl

Note

Integer types (d) read as int (4 bytes), floats (f/g/e) read as double (8 bytes).

Parameters:

• fp -- [in] File pointer to read from (must be opened in binary mode).

• format -- [in] Format string with specifiers (d, f, g, e). Other characters are ignored.

• ... -- [out] Variable pointer arguments matching the format specifiers (int*, double*).

Returns:

The number of items successfully read and assigned. Stops reading on first failure or EOF.

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^2)\) 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).

See also

direct_gaussian_convolution

See also

fft_gaussian_convolution

See also

debug_direct_convolution

Parameters:

• density_grid -- [in] Input density grid array (values corresponding to log_r_grid).

• grid_size -- [in] Number of points in the input grid and density arrays.

• log_r_grid -- [in] Array of logarithmic radial grid coordinates (log10(r)). Must be uniform if FFT is used.

• sigma_log -- [in] Width (standard deviation) of the Gaussian kernel in log10-space.

• result -- [out] Output array (pre-allocated, size grid_size) for the smoothed density field.

void generate_ics_hernquist_anisotropic(double **particles, int npts_initial, gsl_rng *rng, double halo_mass, double scale_radius, gsl_spline **splinemass_out, gsl_interp_accel **enclosedmass_out, gsl_spline **splinePsi_out, gsl_interp_accel **Psiinterp_out, gsl_spline **splinerofPsi_out, gsl_interp_accel **rofPsiinterp_out, gsl_interp **fofEinterp_out, gsl_interp_accel **fofEacc_out, double **radius_out, double **radius_unsorted_out, double **mass_out, double **Psivalues_out, int *num_points_out, int splines_only)

Generates initial conditions for an anisotropic Hernquist profile.

Implements a sampling algorithm for the anisotropic Hernquist distribution function \(f(E, L)\). It uses a combination of inverse transform sampling for radius and direction, and rejection sampling for velocity magnitude.

long long get_available_disk_space(const char *path)

Gets the available disk space for the filesystem containing the given path.

This function uses platform-specific APIs to determine the free space available to the current user on the filesystem where path resides. On Windows, it uses GetDiskFreeSpaceEx. On POSIX-compliant systems (Linux, macOS), it uses statvfs.

Parameters:

path -- [in] A path to a file or directory on the filesystem to check. For Windows, this can be a root directory like "C:\\". For POSIX, any path within the target filesystem, e.g., "data/".

Returns:

long long Available disk space in bytes. Returns -1 on error or if the functionality is not implemented for the current platform.

static inline double get_deterministic_rand(uint64_t seed, int step, int id, int call)

Generates a deterministic, thread-safe uniform random number in [0, 1).

This function hashes the unique coordinates of a simulation event (Seed, Time, ParticleID, CallID) to produce a statistically independent random number. The counter-based approach eliminates shared state: each thread-particle-timestep combination maps to a unique hash input, ensuring reproducibility regardless of thread scheduling or execution order.

Algorithm: Combines the four input coordinates through bit-mixing operations, applies a 64-bit hash finalizer, and converts the result to a uniform double by extracting the high 53 bits and scaling to [0, 1).

Parameters:

• seed -- [in] Simulation seed for reproducibility.

• step -- [in] Current timestep index.

• id -- [in] Persistent particle identifier.

• call -- [in] RNG call sequence number within timestep.

Returns:

double Uniform random value in [0, 1).

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 (e.g., "quadsort_parallel") to user-friendly descriptive names (e.g., "Parallel Quadsort"). This is primarily used for display purposes in command-line output or logs, providing more context than the internal short string identifiers.

Parameters:

sort_alg -- [in] The internal algorithm identifier string.

Returns:

const char* A descriptive name for the algorithm. If no match is found, the input sort_alg string itself is returned as a fallback.

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.

Parameters:

• base_filename -- [in] Original filename.

• with_suffix -- [in] Flag indicating whether to apply the suffix (1=yes, 0=no).

• buffer -- [out] Output buffer for the resulting filename.

• bufsize -- [in] Size of the output buffer.

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.

Note

Returns 0.0 if use_identity_gravity is set to 1.

Note

\(M(r)\) is approximated as (current_rank / npts) * halo_mass_value.

Parameters:

• r -- [in] Radius in kpc.

• current_rank -- [in] Particle rank (0 to npts-1), used for \(M(r)\) approximation.

• npts -- [in] Total number of particles.

• G_value -- [in] Gravitational constant value (e.g., G_CONST).

• halo_mass_value -- [in] Total halo mass (e.g., HALO_MASS).

Returns:

Gravitational acceleration (force per unit mass) in simulation units (kpc/Myr \(^2\)).

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.

See also

gravitational_force

See also

doLeviCivitaLeapfrog

Note

Returns 0.0 if use_identity_gravity is set to 1.

Parameters:

• rho -- [in] Transformed radial coordinate (sqrt(r)). Units: sqrt(kpc).

• current_rank -- [in] Particle rank (0 to npts-1).

• npts -- [in] Total number of particles.

• G_value -- [in] Gravitational constant value (e.g., G_CONST).

• halo_mass_value -- [in] Total halo mass (e.g., HALO_MASS).

Returns:

\(dv/d\tau\) in Levi-Civita coordinates (kpc \(^2\)/Myr \(^2\)).

static void handle_sidm_step(double **particles, int npts, double dt, double current_sim_time, double active_profile_rc, int current_method_display_num, int bootstrap_phase_active, int *current_scatter_counts, int timestep_idx)

Handles the SIDM scattering phase for a single timestep.

Checks if SIDM is enabled and if not in a bootstrap phase that should skip SIDM. If proceeding, it resets particle scatter flags, selects serial or parallel execution based on g_sidm_execution_mode, calls the appropriate core scattering function (perform_sidm_scattering_serial or perform_sidm_scattering_parallel), updates the global total scatter count, and logs debug information if scatters occurred and debugging is enabled. The core scattering functions are responsible for updating the g_particle_scatter_state flags for particles that underwent scattering.

Note

This function modifies the particles array in-place.

Note

It uses global variables: g_enable_sidm_scattering, g_sidm_execution_mode, g_sidm_seed, g_total_sidm_scatters, g_active_halo_mass, g_doDebug, and g_particle_scatter_state.

Parameters:

• particles -- [in,out] The main particle data array: particles[component][current_sorted_index]. Modified in-place with post-scattering velocities/angular momenta.

• npts -- [in] Total number of particles.

• dt -- [in] The simulation timestep (Myr).

• current_sim_time -- [in] The current simulation time at the beginning of this step (Myr).

• active_profile_rc -- [in] The scale radius (kpc) of the currently active profile (NFW, Cored, or Hernquist), passed to sigmatotal.

• current_method_display_num -- [in] The user-facing display number of the current integration method (for logging).

• bootstrap_phase_active -- [in] Flag (0 or 1) indicating if a bootstrap phase (e.g., for Adams-Bashforth) is active. If 1, SIDM scattering is skipped for this step.

double hernquist_potential_wrapper(double r, void *params)

Wrapper function for Hernquist potential calculation.

Parameters:

• r -- Radius (kpc).

• params -- Void pointer to hernquist_potential_params struct.

Returns:

Potential value at radius r.

static inline double hyperg_2F1_safe(double beta, double E_tilde)

Safe wrapper for hypergeometric function with analytical special cases.

Handles boundary cases where GSL has singularities:

• \(\beta=0.5\): \(a=0\), so \({}_2F_1(0,b,c,x) = 1\) (exact)

• \(\beta=-0.5\): \({}_2F_1(2,6,4,x) = (3x^2-10x+10)/(5(x-1)^2)\) (analytical)

• Otherwise: Uses GSL implementation

Parameters:

• beta -- Anisotropy parameter

• E_tilde -- Dimensionless energy (must be < 1 for convergence)

Returns:

Hypergeometric function value

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.

Parallel implementation of insertion sort algorithm optimized for particle sorting.

Uses multiple OpenMP threads to sort sections of particle data in parallel, followed by seam-fixing operations to ensure global ordering. Includes dynamic section calculation and optimized overlap sizing based on chunk characteristics.

Parameters:

• columns -- 2D array of particle data to be sorted

• n -- Number of elements to sort

• columns -- Column-major data array [particle][component]

• n -- Number of particles to sort

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 insertion_sort_sub(double **columns, int start, int end)

Performs insertion sort on a subarray of particle data.

Operates on range [start..end] inclusive within the columns array. Used by 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).

static int isFloat(const char *str)

Checks if a given string represents a valid floating-point number.

Validates if the input string conforms to common floating-point number formats, including an optional leading sign ('+' or '-'), digits, at most one decimal point (if not in exponent part), and an optional exponent part (e.g., "e+10", "E-5"). The function requires at least one digit to be present for a number to be considered valid (e.g., "." or "+." are not valid floats).

Parameters:

str -- [in] The null-terminated string to check.

Returns:

int 1 if the string is a valid float, 0 otherwise.

static int isInteger(const char *str)

Checks if a given string represents a valid integer.

Allows an optional leading '+' or '-' sign. Validates that all subsequent characters in the string are digits. Returns 0 (false) for empty strings, strings containing only a sign, or strings with non-digit characters after the optional sign.

Parameters:

str -- [in] The null-terminated string to check.

Returns:

int 1 if the string is a valid integer, 0 otherwise.

static int load_chosen_particles(const char *filename, int **chosen_out, int *nlowest_out, int **traj_ids_out, int *n_traj_out)

Load chosen particle IDs from a binary file.

This function loads particle IDs that were previously saved for trajectory tracking, restoring the particle selection state from a prior simulation run. It reads two sets of particle IDs: one for particles with the lowest angular momentum, and another for general trajectory tracking by original particle ID.

The expected binary file format is:

• Header (8 bytes):

  • int32: num_lowestl - Number of particles with lowest angular momentum

  • int32: num_trajectories - Number of particles for general trajectory tracking

• Data section:

  • int32[num_lowestl]: Particle IDs selected by lowest angular momentum

  • int32[num_trajectories]: Particle IDs for general trajectory tracking

For backward compatibility, the function also supports files with only a single integer header followed by lowest_l particle IDs. This allows older simulation outputs to be read correctly.

The loaded particle IDs are used to maintain continuity in trajectory tracking across simulation restarts and extensions. The lowest_l particles continue to be tracked in lowest_l_trajectories.dat, while the trajectory particles continue in trajectories.dat and single_trajectory.dat.

Note

If traj_ids_out is NULL, trajectory IDs are read but not returned to the caller. Memory is allocated for the output arrays; caller is responsible for freeing them.

Parameters:

• filename -- [in] Path to the input binary file containing saved particle IDs

• chosen_out -- [out] Pointer to store allocated array of lowest_l particle IDs

• nlowest_out -- [out] Number of lowest_l particles loaded

• traj_ids_out -- [out] Pointer to store allocated array of trajectory particle IDs (can be NULL)

• n_traj_out -- [out] Number of trajectory particles loaded (can be NULL)

Returns:

1 on success, 0 on failure (file not found or read error)

static int load_existing_energy_diagnostics(const char *filename, int max_snapshots)

Writes collected debug energy data to a file.

Writes the time evolution of theoretical vs. dynamic energy for the tracked particle (DEBUG_PARTICLE_ID) to data/debug_energy_compare.dat. The output includes kinetic and potential energy components. This function is called at the end of the simulation if g_doDebug is enabled.

Loads existing energy diagnostic data from a file for restart/extend operations.

See also

store_debug_approxE

See also

store_debug_dynE_components

Reads energy data from an existing total_energy_vs_time file to populate the global energy arrays. This ensures continuity when restarting simulations.

Note

File writing is suppressed if skip_file_writes is non-zero.

Parameters:

• filename -- [in] Path to the existing energy file

• max_snapshots -- [in] Maximum number of snapshots to load

Returns:

Number of snapshots successfully loaded

static int load_from_double_buffer(const char *tag, int snapshot_index, int total_snapshots, double **particles, int npts, int *particle_ids)

Attempts to load a snapshot from the double-precision buffer files.

Tries to load the specified snapshot from the double buffer files. Returns 1 if successful, 0 if snapshot not in buffer or files missing. On success, loads particle state with full double precision (no conversion).

Buffer stores snapshots in circular order. For 81-snapshot run, buffer contains snapshots [77, 78, 79, 80]. Snapshot 79 is located at buffer position 2.

Parameters:

• tag -- Simulation tag for filename

• snapshot_index -- Snapshot number to load (e.g., 79)

• total_snapshots -- Total snapshots written so far (e.g., 81)

• particles -- Destination particle arrays [R, Vrad, L, ID, \(\mu\), cos_phi, sin_phi]

• npts -- Number of particles

• particle_ids -- Destination for particle IDs (unused in this function).

Returns:

1 if loaded successfully, 0 if snapshot not available in buffer

static void load_particles_from_restart(const char *filename, int snapshot_index, double **particles, int npts, int block_size, int *inverse_map)

Load particle state from all_particle_data file for restart/extend operations.

Reads a specific snapshot from an all_particle_data file and converts it to the internal particles array format used for initial conditions. This function is used by --sim-restart and --sim-extend to load the final particle state from a completed simulation run.

Parameters:

• filename -- Path to the all_particle_data file

• snapshot_index -- Which snapshot to load (0-based)

• particles -- Output array [5][npts] in standard IC format

• npts -- Number of particles

• block_size -- Block size for I/O (typically 100)

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

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

See also

g_enable_logging

Note

Creates the "log" directory if it does not exist.

Warning

Prints an error to stderr if the log file cannot be opened. Logging only occurs if the global g_enable_logging flag is set.

Parameters:

• level -- [in] Severity level (e.g., "INFO", "WARNING", "ERROR").

• format -- [in] Printf-style format string.

• ... -- [in] Variable arguments for the format string.

threevector make_threevector(double x, double y, double z)

Constructs a three-dimensional vector from its Cartesian components.

This utility function initializes a threevector structure with the provided x, y, and z components. It serves as a convenient constructor.

Parameters:

• x -- [in] The x-component of the vector.

• y -- [in] The y-component of the vector.

• z -- [in] The z-component of the vector.

Returns:

threevector An initialized threevector structure.

double mass_integrand_hernquist(double r, void *p)

GSL-compatible wrapper for the Hernquist density integrand: \(4\pi r^2 \rho(r)\).

Parameters:

• r -- Radius (kpc).

• p -- Void pointer to hernquist_params struct.

Returns:

Integrand value for mass calculation.

double massintegrand(double r, void *params)

GSL integrand \(r^2 \rho_{shape}(r)\) for Cored Plummer-like profile mass calculation.

Computes the term \(r^2 \rho_{shape}(r)\) for the Cored Plummer-like density profile, where \(\rho_{shape}(r) = (1 + (r/r_c)^2)^{-3}\). This integrand is used in GSL numerical integration routines (e.g., gsl_integration_qag) to calculate the normalization factor or the enclosed mass \(M(<r) = 4\pi \int_0^r r'^2 \rho_{physical}(r') dr'\). It directly uses the RC macro for the scale radius.

Parameters:

• r -- [in] Radial coordinate \(r\) (kpc).

• params -- [in] Void pointer to parameters (unused).

Returns:

double The value of the mass integrand \(r^2 \rho_{shape}(r)\) at radius r.

double massintegrand_hernquist(double r, void *params)

GSL integrand \(r^2 \rho(r)\) for Hernquist profile mass calculation.

Computes \(r^2 \rho(r)\) where \(\rho(r) = M_{\text{total}} a / [2\pi r (r+a)^3]\) is the Hernquist density profile.

Parameters:

• r -- [in] Radial coordinate (kpc).

• params -- [in] Void pointer to a double array: [a_scale, M_total, normalization_factor].

Returns:

double The value of \(r^2 \rho(r)\) at radius r.

double massintegrand_profile_nfwcutoff(double r, void *params)

GSL integrand \(r^2 \rho(r)\) for NFW-like profile mass calculation.

Computes \(r^2 \rho(r)\) where \(\rho(r)\) is an NFW-like density profile with an inner softening term and an outer power-law cutoff. The density \(\rho(r) = p[2] \times [ (r_s + \epsilon)(1+r_s)^2 (1 + (r_s/C)^N) ]^{-1}\), where \(r_s = r/p[0]\), \(\epsilon=0.01\), \(C=p[3]\), \(N=10.0\). This integrand is used in GSL routines to calculate the normalization factor or enclosed mass for the NFW-like profile.

Parameters:

• r -- [in] Radial coordinate \(r\) (kpc).

• params -- [in] Void pointer to a double array p of size 4:

  • p[0] (rc_param): Scale radius (kpc).

  • p[1] (halo_mass): Target total halo mass ( \(M_{\odot}\)) - used to derive nt_nfw_scaler.

  • p[2] (nt_nfw_scaler): Density normalization constant \(nt_{NFW}\).

  • p[3] (falloff_C_param): Falloff transition factor \(C\).

Returns:

double The value of the mass integrand \(r^2 \rho(r)\). Returns 0.0 if rc_param (p[0]) is non-positive.

static inline uint64_t mix64(uint64_t k)

High-performance 64-bit mixing function.

Applies three rounds of reversible bit-mixing operations (XOR-shift and multiply) to decorrelate input bits. This mixing pattern provides excellent avalanche properties: single-bit changes in the input propagate to affect approximately half of the output bits, ensuring uniform distribution.

Parameters:

k -- [in] Input 64-bit integer key.

Returns:

uint64_t The mixed 64-bit hash value.

double nfw_potential_wrapper(double r, void *params)

Wrapper function for NFW potential calculation.

Parameters:

• r -- Radius (kpc).

• params -- Void pointer to nfw_potential_params struct.

Returns:

Approximate potential value at radius r.

double om_mu_integrand(double mu, void *params)

Integrand for OM \(\mu\)-integral in df_fixed_radius.

static void parallel_radix_sort(double **columns, int n)

CPU-based parallel radix sort optimized for particle data.

High-performance radix sort that operates directly on double-precision radius values by converting them to sortable integer representations. Automatically chooses between serial and parallel implementations based on array size. Typically outperforms GPU sorting for arrays < 50M elements due to lower overhead.

Parameters:

• columns -- 2D array of particle data to be sorted [particle][component]

• n -- Number of particles to sort

int parse_nsphere_filename(const char *filename, int *N, int *Ntimes, double *tfinal)

Parse NSphere filename to extract simulation parameters.

Attempts to extract N, Ntimes, and tfinal from standard NSphere filename format: prefix_tag_N_Ntimes_tfinal.dat (e.g., beta025_1000_1250000_100.dat)

Parameters:

• filename -- The filename to parse

• N -- Pointer to store extracted N value

• Ntimes -- Pointer to store extracted Ntimes value

• tfinal -- Pointer to store extracted tfinal value

Returns:

1 if successful, 0 if unable to parse

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

Parses sub-arguments for the --save command-line option.

This function is called when the --save option is encountered during command-line argument parsing. It reads subsequent arguments (until another option starting with '-' is found, or arguments end) which specify the level or type of data to save. It then sets the corresponding global data output flags (g_doDebug, g_doDynPsi, g_doDynRank, g_doAllParticleData) based on the highest priority valid sub-argument encountered. Valid sub-arguments and their priority (lowest to highest):

• "raw-data": Enables g_doAllParticleData.

• "psi-snaps": Enables g_doAllParticleData, g_doDynPsi.

• "full-snaps": Enables g_doAllParticleData, g_doDynPsi, g_doDynRank.

• "all" or "debug-energy": Enables all flags (g_doDebug, g_doDynPsi, g_doDynRank, g_doAllParticleData).

Note

Exits the program with an error message if an unknown sub-argument to --save is found.

Parameters:

• argc -- [in] The total argument count from main().

• argv -- [in] The argument array from main().

• pIndex -- [in,out] Pointer to the current index in argv. On input, it points to the --save option. On output, it is updated to point to the last sub-argument consumed by this function.

static inline void perform_scatter_update(double **particles, int i, int j, uint64_t seed, int step, int pid)

void perform_sidm_scattering_parallel_graphcolor (double **__restrict__ particles, int npts, double dt, int timestep_idx, uint64_t sidm_seed, long long *Nscatter_total_step, double halo_mass_for_sidm, double rc_for_sidm, int *current_scatter_counts)

Helper function to perform scatter update for two particles.

Updates velocities and angular components for both particles after a scattering event. The graph coloring algorithm uses this function to directly update particle states, avoiding intermediate buffering.

Performs SIDM scattering calculations for one timestep using graph coloring parallel algorithm.

Implements a parallel SIDM scattering algorithm using graph coloring to eliminate race conditions. Particles are divided into 11 color groups where particles of the same color are separated by more than the interaction range (>10 positions), ensuring no conflicts. Each color is processed sequentially, but particles within a color are processed in parallel with direct state updates and no synchronization overhead.

Azimuthal Angle Treatment: Uses persistent azimuthal angle storage (cos(φ), sin(φ)) in particles[5] and particles[6]. These values are initialized randomly at simulation start and updated only during scattering events. Between scatters, φ remains constant (no orbital precession tracking). This approximation is physically valid because in high-scatter regimes φ evolves minimally between events, while in low-scatter regimes stochastic partner selection and orbital dynamics erase correlations from prior scatters. This approach maintains accuracy while avoiding costly per-timestep random number generation.

Graph Coloring Algorithm:

1. Divides particles into 11 colors (modulo arithmetic: color = i % 11).

2. Processes each color sequentially to prevent neighbor conflicts.

3. Within each color, processes particles in parallel with per-thread RNGs.

4. For each particle: a. Loads persistent φ orientation from stored cos(φ), sin(φ). b. Evaluates scattering probability with neighbors using 3D relative velocities. c. If scatter occurs: performs isotropic scattering, updates velocities and φ storage. d. Updates g_particle_scatter_state for AB3 integrator compatibility.

This design achieves high parallel efficiency without race conditions or atomic operations.

Parameters:

• particles -- Particle data arrays

• i -- First particle index

• j -- Second particle index

• rng -- Random number generator for scattering angles

• particles -- [in,out] Main particle data array: particles[component][current_sorted_index]. Modified in-place with post-scattering velocities/angular momenta.

• npts -- [in] Total number of particles.

• dt -- [in] Simulation timestep (Myr), used in probability calculation.

• current_time -- [in] Current simulation time (Myr). Marked unused but available for future use.

• rng_per_thread_list -- [in] Array of GSL RNG states, one for each OpenMP thread.

• num_threads_for_rng -- [in] The number of allocated RNGs in rng_per_thread_list (equals max threads).

• Nscatter_total_step -- [out] Pointer to a long long to accumulate the total number of scatter events that occurred in this timestep.

• halo_mass_for_sidm -- [in] Total halo mass ( \(M_{\odot}\)) for the active profile, passed to sigmatotal.

• rc_for_sidm -- [in] Scale radius (kpc) for the active profile, passed to sigmatotal.

• current_scatter_counts -- [in,out] Per-particle scatter count array for tracking (can be NULL).

void perform_sidm_scattering_parallel_graphcolor(double **particles, int npts, double dt, int timestep_idx, uint64_t sidm_seed, long long *Nscatter_total_step, double halo_mass_for_sidm, double rc_for_sidm, int *current_scatter_counts)

void perform_sidm_scattering_serial(double **particles, int npts, double dt, int timestep_idx, uint64_t sidm_seed, long long *Nscatter_total_step, double halo_mass_for_sidm, double rc_for_sidm, int *current_scatter_counts)

Executes self-interacting dark matter (SIDM) scattering for one simulation timestep (Serial version).

Implements a serial SIDM scattering algorithm with persistent azimuthal angle tracking. For each particle i (the primary scatterer), it considers up to nscat (typically 10) subsequent particles in the array as potential scattering partners. The particles array is sorted by radius, making these neighbors spatially proximate.

Azimuthal Angle Treatment: The azimuthal angle φ for each particle's perpendicular velocity is stored persistently as cos(φ) and sin(φ) in particles[5] and particles[6]. These values are initialized randomly at simulation start and updated only during scattering events. Between scatters, φ remains constant (no orbital precession tracking). This approximation is valid because:

• In high-scatter regimes, φ evolves minimally between frequent scattering events.

• In low-scatter regimes, stochastic partner selection and orbital dynamics erase correlations from prior scattering history. This approach maintains physical accuracy while avoiding costly random number generation at every timestep.

Algorithm Steps:

1. Loads persistent φ orientation (cos(φ), sin(φ)) from particles[5] and particles[6].

2. For each of the nscat potential partners m: a. Loads partner's φ orientation and calculates relative azimuthal separation. b. Computes 3D relative velocity \(v_{rel}\) using Law of Cosines. c. Calculates interaction rate \(\Gamma_m = \sigma(v_{rel}) v_{rel}\). d. Accumulates total interaction rate \(\Gamma_{tot} = \sum \Gamma_m\).

3. Estimates local particle number density using shell volume approximation.

4. Calculates scattering probability: \(P_{scatter} = 1 - \exp(-\Gamma_{tot} \times dt / (2V_{shell}))\).

5. If scatter occurs (stochastic decision): a. Selects partner weighted by individual \(\Gamma_m\) contributions. b. Performs isotropic scattering in center-of-mass frame. c. Updates both particles' velocities (radial and perpendicular components). d. Updates persistent φ storage (cos(φ), sin(φ)) for both particles. e. Sets g_particle_scatter_state flags for AB3 integrator history management.

Parameters:

• particles -- [in,out] Main particle data array: particles[component][current_sorted_index]. Modified in-place with post-scattering velocities/angular momenta.

• npts -- [in] Total number of simulation particles.

• dt -- [in] Integration timestep (Myr).

• current_time -- [in] Current simulation time (Myr). Marked unused but available.

• rng -- [in] GSL random number generator instance for all stochastic processes.

• Nscatter_total_step -- [out] Pointer to a long long to accumulate total scattering events this timestep.

• halo_mass_for_sidm -- [in] Total halo mass ( \(M_{\odot}\)) for the active profile, passed to sigmatotal.

• rc_for_sidm -- [in] Scale radius (kpc) for the active profile, passed to sigmatotal.

• current_scatter_counts -- [in,out] Per-particle scatter count array for tracking (can be NULL).

static inline double potential_hernquist(double r, double M, double a)

Analytical potential for the Hernquist profile.

Parameters:

• r -- Radius (kpc).

• M -- Total halo mass ( \(M_{\odot}\)).

• a -- Scale radius \(a\) of the Hernquist profile (kpc).

Returns:

Gravitational potential in code units of (kpc/Myr) \(^2\).

static void printUsage(const char *prog)

Displays detailed usage information for command-line arguments.

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.

Note

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

Parameters:

prog -- [in] The program name to display in the usage message (typically argv[0]).

int prompt_yes_no(const char *prompt)

Prompts the user with a yes/no question and reads their response from stdin.

Displays the given prompt string followed by "[y/N]: ". Reads a line of input from the user.

• Returns 1 (yes) if the first character of the input is 'y' or 'Y'.

• Returns 0 (no) if the first character is 'n', 'N', or if the input is an empty line (user just pressed Enter, defaulting to No).

• If any other input is received, the prompt is repeated. Handles potential EOF or read errors by defaulting to No.

Parameters:

prompt -- [in] The question/prompt message to display to the user.

Returns:

int 1 if the user confirms (yes), 0 otherwise (no/default).

double Psiintegrand(double rp, void *params)

Integrand for calculating gravitational potential.

Computes the integrand required for the integral term in the potential calculation, \(\Psi(r)\). This function acts as a dispatcher, calling the appropriate profile-specific mass integrand via a function pointer provided in the params struct.

See also

Psiintegrand_params

See also

massintegrand

See also

massintegrand_profile_nfwcutoff

Parameters:

• rp -- The radial integration variable \(r'\) (kpc).

• params -- Void pointer to a Psiintegrand_params struct, which contains the function pointer to the profile's mass integrand.

Returns:

The value of the potential integrand at rp.

void quadsort_parallel_sort(double **columns, int n)

Parallel quadsort implementation using chunk-based approach with overlap.

Uses quadsort algorithm for both initial chunk sorting and seam fixing, providing better performance than insertion sort for large datasets while maintaining sort stability.

Parallel implementation of quadsort algorithm optimized for particle sorting.

Uses multiple OpenMP threads to sort sections of particle data in parallel, followed by seam-fixing operations to ensure global ordering. Includes dynamic section calculation and optimized overlap sizing based on chunk characteristics.

Parameters:

• columns -- 2D array of particle data to be sorted

• n -- Number of elements to sort

• columns -- Column-major data array [particle][component]

• n -- Number of 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 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").

See also

write_initial_conditions

Note

Expects 7 variables per particle. See write_initial_conditions for file format details.

Warning

Prints an error to stderr if the file cannot be opened, if the particle count does not match npts (returns early), or if a read error occurs.

Parameters:

• particles -- [out] 2D array to store loaded particle properties [component][particle_index]. Must be pre-allocated with dimensions [7][npts].

• npts -- [in] Expected number of particles to read.

• filename -- [in] Path to the input binary file.

static void reassign_orig_ids_with_rank(double *orig_ids, int n)

Remaps particle IDs to zero-based rank order within the current particle set.

Typically used after tidal stripping where a subset of particles remains. The orig_ids array contains potentially non-contiguous IDs of the n remaining particles. Sorts these IDs and replaces each with its rank (0 to n-1) within the sorted sequence. Transforms arbitrary ID values into a compact, contiguous sequence of rank IDs. The input orig_ids array (which is particles[3] in main) is modified in-place.

Note

Uses qsort and a temporary array for sorting. Exits via exit(1) if memory allocation for the temporary array fails.

Parameters:

• orig_ids -- [in,out] Pointer to an array of particle IDs (stored as doubles). These are modified in-place to become rank IDs.

• n -- [in] The number of elements in the orig_ids array (i.e., npts after stripping).

static void retrieve_all_particle_phi_snapshot(const char *filename, int snap, int npts, int block_size, float *phi_out)

Retrieves particle phi angles for a specific snapshot from the all_particle_phi.dat binary file.

This function reads the phi angle data for all npts particles corresponding to a single snapshot number (snap) from the specified binary file. The file is expected to be in step-major order, where each record per particle consists of a single float (phi). It calculates the correct file offset to seek to the desired snapshot. To ensure thread safety when called in parallel (e.g., during post-processing of snapshots), file I/O (seeking and reading) is performed within an OpenMP critical section named file_access_phi. Temporary local buffers are used for reading, and data is then copied to the caller-provided output array.

Note

Exits via CLEAN_EXIT(1) on memory allocation failure, file open failure, fseek failure, or unexpected EOF.

Parameters:

• filename -- [in] Path to the binary phi data file (e.g., "data/all_particle_phi<suffix>.dat").

• snap -- [in] The snapshot number (0-indexed, corresponding to write events) to retrieve.

• npts -- [in] Number of particles per snapshot.

• block_size -- [in] Number of snapshots written per block to file.

• phi_out -- [out] Pointer to an array (size npts) to store the retrieved phi angle values.

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 the all_particle_data.dat binary file.

This function reads the data (rank, radius, radial velocity, angular momentum) for all npts particles corresponding to a single snapshot number (snap) from the specified binary file. The file is expected to be in step-major order, where each record per particle consists of an int (rank) and three floats (R, Vrad, L). It calculates the correct file offset to seek to the desired snapshot. To ensure thread safety when called in parallel (e.g., during post-processing of snapshots), file I/O (seeking and reading) is performed within an OpenMP critical section named file_access. Temporary local buffers are used for reading, and data is then copied to the caller-provided output arrays.

Note

Exits via CLEAN_EXIT(1) on memory allocation failure, file open failure, fseek failure, or unexpected EOF.

Parameters:

• filename -- [in] Path to the binary data file (e.g., "data/all_particle_data<suffix>.dat").

• snap -- [in] The snapshot number (0-indexed, corresponding to write events) to retrieve.

• npts -- [in] Number of particles per snapshot.

• block_size -- [in] Number of snapshots written per block by append_all_particle_data_chunk_to_file. Determines seek offset calculation.

• L_out -- [out] Pointer to an array (size npts) to store the retrieved angular momentum values.

• Rank_out -- [out] Pointer to an array (size npts) to store the retrieved particle ranks.

• R_out -- [out] Pointer to an array (size npts) to store the retrieved radial positions.

• Vrad_out -- [out] Pointer to an array (size npts) to store the retrieved radial velocities.

static void save_chosen_particles(const char *filename, int *chosen, int nlowest, int *traj_ids, int n_traj)

Save chosen particle IDs to a binary file.

This function saves particle IDs for multiple trajectory tracking systems to ensure continuity across restart and extend operations. It writes two sets of particle IDs: one for particles selected based on lowest angular momentum, and another for particles tracked by their original ID for general trajectory analysis.

The binary file format consists of:

• Header (8 bytes):

  • int32: num_lowestl - Number of particles with lowest angular momentum

  • int32: num_trajectories - Number of particles for general trajectory tracking

• Data section:

  • int32[num_lowestl]: Particle IDs selected by lowest angular momentum

  • int32[num_trajectories]: Particle IDs for general trajectory tracking

The lowest_l particles are used by lowest_l_trajectories.dat which tracks radius, energy, and angular momentum for these specific particles over time. The trajectory particles are used by trajectories.dat which tracks radius, velocity, and mu for the specified particles by their original ID. The file single_trajectory.dat uses only the first particle ID from the trajectory list.

This file is essential for maintaining consistent particle tracking when simulations are restarted or extended, ensuring that the same particles continue to be monitored.

Note

Exits via CLEAN_EXIT(1) if file creation fails.

Parameters:

• filename -- [in] Path to the output binary file (e.g., "data/chosen_particles_<suffix>.dat")

• chosen -- [in] Array of particle IDs that have the lowest angular momentum values

• nlowest -- [in] Number of lowest angular momentum particles to save

• traj_ids -- [in] Array of particle IDs to track in trajectories.dat

• n_traj -- [in] Number of trajectory particles to save

double sigmatotal(double vrel, int npts, double halo_mass_for_calc, double rc_for_calc)

Calculates the self-interacting dark matter (SIDM) scattering cross-section.

Implements a velocity-independent cross-section model where the opacity \(\sigma/m = \kappa\) is constant. The function uses the global variable g_sidm_kappa for the opacity \(\kappa\). The total cross-section \(\sigma\) scales with the individual particle mass, \(m_{particle}\), derived from the total halo_mass_for_calc and the number of particles npts. The function includes unit conversions to return the cross-section in simulation units (kpc \(^2\)).

Unit Conversion Detail: \(\kappa\) (cm²/g) \(\times m_{particle}\) ( \(M_{\odot}\)) \(\rightarrow \sigma\) (kpc²) The conversion factor is \(2.089 \times 10^{-10} \text{ (kpc}^2 \text{ } M_{\odot}^{-1}) / (\text{cm}^2 \text{ g}^{-1})\).

Parameters:

• vrel -- [in] Relative velocity between particles (kpc/Myr). This parameter is unused in the velocity-independent model.

• npts -- [in] Total number of simulation particles, used for calculating \(m_{particle}\).

• halo_mass_for_calc -- [in] Total halo mass ( \(M_{\odot}\)) for the active profile, used for \(m_{particle}\).

• rc_for_calc -- [in] Scale radius (kpc) of the active profile. This parameter is unused in the velocity-independent model.

Returns:

double Total scattering cross-section \(\sigma\) in kpc \(^2\). Returns 0.0 if npts or the calculated particle_mass_Msun is non-positive.

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 radial position (rad).

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.

Uses the standard library qsort function with compare_partdata_by_rad as the comparison function. Includes basic safety checks for NULL array or non-positive npts. The sort is performed in-place.

Parameters:

• array -- Array of particle data structures to be sorted

• npts -- Number of particles in the array

• array -- [in,out] Array of PartData structures to be sorted.

• npts -- [in] Number of elements in the array.

void sort_particles(double **particles, int npts)

Convenience wrapper function for sorting particles with the default algorithm.

Calls sort_particles_with_alg using the default sorting algorithm specified in g_defaultSortAlg.

Parameters:

• particles -- [in,out] 2D array of particle data to be sorted [component][particle].

• npts -- [in] Number of particles to sort.

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

Main function for sorting particle data using a specified algorithm.

Transposes the data format from particles[component][particle] to columns[particle][component], applies the specified sorting algorithm, and then transposes back to the original format.

Sorts particle data using the specified sorting algorithm.

Performs a three-phase particle sorting operation:

1. Data transposition to a persistent, column-major buffer.

2. Application of the selected sorting algorithm on the buffer.

3. Reverse transposition of the sorted data back to the original array.

The function uses a persistent global buffer (g_sort_columns_buffer) to minimize memory allocation/deallocation overhead across repeated sort calls.

Parameters:

• particles -- 2D array of particle data to be sorted [component][particle]

• npts -- Number of particles to sort

• sortAlg -- String identifier of the sorting algorithm to use ("quadsort", "quadsort_parallel", or default to "insertion_parallel")

• particles -- 2D array of particle data to be sorted [component][particle]

• npts -- Number of particles to sort

• sortAlg -- Sorting algorithm to use ("quadsort", "quadsort_parallel", "insertion", or "insertion_parallel")

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

Sorts radius and potential arrays in tandem, maintaining their correspondence.

This function takes an array of radial coordinates (rrA_spline) and an array of corresponding potential values (psiAarr_spline). It sorts rrA_spline in ascending order and applies the identical swaps to psiAarr_spline, ensuring that psiAarr_spline[i] still corresponds to rrA_spline[i] after sorting. This is crucial for creating GSL splines where the x-array must be strictly monotonic and the y-array must maintain its pairing with the x-values. The arrays are assumed to have npts elements.

Parameters:

• rrA_spline -- [in,out] Array of radial coordinates to be sorted. Modified in-place.

• psiAarr_spline -- [in,out] Array of corresponding Psi values. Modified in-place in tandem with rrA_spline.

• npts -- The number of points in the arrays (arrays are of size npts).

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 store_debug_approxE(int snapIndex, double E_value, double time_val)

Records the theoretical model energy for a debug snapshot.

See also

store_debug_dynE_components

See also

finalize_debug_energy_output

Note

Updates dbg_count if snapIndex extends the recorded range.

Parameters:

• snapIndex -- [in] Index of the snapshot (0 to DEBUG_MAX_STEPS - 1).

• E_value -- [in] Theoretical energy value (per unit mass).

• time_val -- [in] Simulation time (Myr).

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.

See also

store_debug_approxE

See also

finalize_debug_energy_output

Note

Updates dbg_count if snapIndex extends the recorded range.

Parameters:

• snapIndex -- [in] Index of the snapshot (0 to DEBUG_MAX_STEPS - 1).

• totalE -- [in] Total energy (KE + PE) per unit mass.

• kinE -- [in] Kinetic energy component (per unit mass).

• potE -- [in] Potential energy component (per unit mass).

• time_val -- [in] Simulation time (Myr).

• radius_val -- [in] Particle radius (kpc).

static int truncate_files_to_snapshot(int last_snapshot, int npts, int dtwrite, int num_traj_particles, int nlowest, const char *file_suffix)

Truncate all output files to the last common complete snapshot.

Ensures all files are consistent for restart by removing partial data.

Parameters:

• last_snapshot -- The last complete snapshot to keep (0-based)

• npts -- Number of particles

• dtwrite -- Timestep interval for writes

• num_traj_particles -- Number of trajectory particles

• nlowest -- Number of lowest-L particles

• file_suffix -- The file suffix string

Returns:

0 on success, -1 on error

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

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

Writes particle initial conditions to a binary file.

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.

See also

read_initial_conditions

Note

The binary file format is: npts (int32), followed by npts records, each consisting of 7 double values: (radius, velocity, ang. mom., orig. index, mu, cos_phi, sin_phi). cos_phi and sin_phi are the cosine and sine of the azimuthal angle used for SIDM scattering.

Warning

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

Parameters:

• particles -- [in] 2D array containing particle properties [component][particle_index]. Expected components: 0=rad, 1=vel, 2=angmom, 3=orig_idx, 4=mu.

• npts -- [in] Number of particles to write.

• filename -- [in] Path to the output binary file.

Variables

static int *chosen = NULL

Array of indices (original IDs) for selected low-L particles.

static double dbg_approxE[DEBUG_MAX_STEPS]

Arrays for tracking energy components through simulation for debugging.

Theoretical model energy (per unit mass).

static int dbg_count = 0

Tracks the number of recorded debug snapshots.

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_radius[DEBUG_MAX_STEPS]

Particle radius at each snapshot (kpc).

static double dbg_time[DEBUG_MAX_STEPS]

Simulation time at each snapshot (Myr).

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 doReadInit = 0

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

static int doWriteInit = 0

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

static double g_active_halo_mass = HALO_MASS

Active halo mass for N-body force calculations.

static double g_anisotropy_beta = 0.0

Anisotropy parameter \(\beta\) for Hernquist profile.

static int g_anisotropy_beta_provided = 0

Flag: 1 if the --aniso-beta option is provided.

static int g_attempt_load_seeds = 0

Flag: 1 to attempt loading seeds from files if not provided via command line.

static int g_benchmark_count = 0

Number of times benchmarking was performed.

static double g_cored_profile_halo_mass = HALO_MASS

Total halo mass for the Cored profile ( \(M_{\odot}\)). Populated from g_halo_mass_param.

static double g_cored_profile_rc = RC

Scale radius \(r_c\) for the Cored profile (kpc). Populated from g_scale_radius_param.

static double g_cored_profile_rmax_factor = CUTOFF_FACTOR_CORED_DEFAULT

Cutoff radius factor for the Cored profile. Populated from g_cutoff_factor_param.

static char g_current_best_sort[32] = "insertion_parallel"

Tracks the best-performing sorting algorithm for adaptive mode.

static int *g_current_timestep_scatter_counts = NULL

Tracks scatter counts per particle for the timestep block being processed. The counts are reset after each block write.

static double g_cutoff_factor_param = CUTOFF_FACTOR_CORED_DEFAULT

Generalized rmax factor, defaults to Cored's default.

static int g_cutoff_factor_param_provided = 0

Flag: 1 if the --cutoff-factor option is provided.

static int g_dbl_buf_count = 0

Number of valid snapshots in buffer (range: 0-4, saturates at 4).

static int g_dbl_buf_current_slot = 0

Next write slot in rolling 4-snapshot circular buffer (range: 0-3, wraps via modulo).

static double *g_dbl_buf_L = NULL

Buffer for particle angular momenta (last 4 snapshots, npts particles each).

static int g_dbl_buf_npts = 0

Number of particles, used for buffer indexing.

static double *g_dbl_buf_phi = NULL

Buffer for particle phi angles (last 4 snapshots, npts particles each).

static double *g_dbl_buf_R = NULL

Rolling buffer system for preserving last 4 snapshots in full double precision.

Maintains a circular buffer of the most recent snapshots to avoid precision loss during restart/extend operations. The buffer stores:

• Snapshots 0, then [0,1], then [0,1,2], then [0,1,2,3]

• After 4 snapshots: rolling last 4 (e.g., [77,78,79,80] for 81-snapshot run)

Files written: double_buffer_all_particle_data_<tag>.dat (28 bytes/particle/snap) double_buffer_all_particle_phi_<tag>.dat (8 bytes/particle/snap)

Buffer size: (28 + 8) * npts * 4 = 144 * npts bytes in memory Buffer for particle radii (last 4 snapshots, npts particles each).

static int *g_dbl_buf_Rank = NULL

Buffer for particle ranks (last 4 snapshots, npts particles each).

static double *g_dbl_buf_Vrad = NULL

Buffer for particle radial velocities (last 4 snapshots, npts particles each).

static int g_default_max_threads = 0

Default max threads (all cores) for non-SIDM operations.

static const char *g_defaultSortAlg = "quadsort_parallel"

< Default sorting algorithm identifier string. Set based on command-line options.

static drhodr_func_t g_density_derivative_func = NULL

Function pointer for dynamically selected density derivative (Cored/Hernquist profiles).

Points to the appropriate \(d\rho/dr\) implementation selected at runtime during IC generation. Enables profile-specific calculations without conditional branching in performance-critical loops.

Selection logic:

• If OM model active: Points to augmented density derivative (density_derivative_om_cored)

• Otherwise: Points to standard isotropic derivative (drhodr)

Note

Must be set before any integration begins.

static drhodr_nfw_func_t g_density_derivative_nfw_func = NULL

Function pointer for dynamically selected NFW density derivative.

Points to the appropriate NFW-specific \(d\rho/dr\) implementation selected at runtime during IC generation. NFW derivatives require additional parameters (rc, nt_nfw, falloff_C) passed through the function signature.

Selection logic:

• If OM model active: Points to augmented NFW derivative (density_derivative_om_nfw)

• Otherwise: Points to standard NFW derivative (drhodr_profile_nfwcutoff)

Note

Must be set before any integration begins.

int g_doAllParticleData = 1

Save complete particle evolution history (default: on).

int g_doDebug = 1

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 (default: on).

int g_doDynPsi = 1

Enable dynamic potential recalculation (default: on).

int g_doDynRank = 1

Enable dynamic rank calculation per step (default: on).

int g_doRestart = 0

Enable simulation restart from checkpoint.

int g_doRestartForce = 0

Force regeneration of all snapshots when restarting.

int g_doSimExtend = 0

Enable simulation extension mode (--sim-extend).

int g_doSimRestart = 0

Enable simulation restart detection mode (--sim-restart).

int g_doSimRestartCheckOnly = 0

Check-only mode for --sim-restart (do not actually restart).

static double **g_E_buf = NULL

Buffer for total relative energy E_rel [particle_idx][timestep_in_buffer].

static double *g_E_init_vals = NULL

Stores initial energy E_rel for tracked particles [particle_idx].

int g_enable_logging = 0

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

int g_enable_sidm_scattering = 0

Enable SIDM scattering physics (0=no, 1=yes). Default is OFF.

static FILE *g_energy_file = NULL

File handle for energy_and_angular_momentum_vs_time.dat.

static int g_energy_snapshots_loaded = 0

Number of energy snapshots loaded from existing file in restart mode.

char *g_extend_file_source = NULL

Source file to extend from (--extend-file).

static double g_falloff_factor_param = FALLOFF_FACTOR_NFW_DEFAULT

Generalized falloff factor, defaults to NFW's default.

static int g_falloff_factor_param_provided = 0

Flag: 1 if the --falloff-factor option is provided.

static char g_file_suffix[256] = ""

Global file suffix string.

static double g_halo_mass_param = HALO_MASS

Generalized halo mass ( \(M_{\odot}\)), defaults to HALO_MASS macro.

static int g_halo_mass_param_provided = 0

Flag: 1 if the --halo-mass option is provided.

static int g_hybrid_p_cores = 0

Number of P-cores on hybrid CPUs (0 if not hybrid).

static unsigned long int g_initial_cond_seed = 0

Seed used for generating initial conditions.

static const char *g_initial_cond_seed_filename_base = "data/last_initial_seed"

Base name for IC seed file.

static int g_initial_cond_seed_provided = 0

Flag: 1 if the --init-cond-seed option is provided.

static int g_insertion_wins = 0

Number of times insertion sort was faster.

static double **g_L_buf = NULL

Buffer for angular momentum [particle_idx][timestep_in_buffer].

static double *g_L_init_vals = NULL

Stores initial angular momentum \(L\) for tracked particles [particle_idx].

static double g_l_target_value = -1.0

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 \(L_{compare}\). Target \(L\) value for --lvals-target mode. A negative value indicates the option was not set.

static double **g_lowestL_E_buf = NULL

Buffer for energy of low-L particles [particle_idx][timestep_in_buffer].

static FILE *g_lowestL_file = NULL

File handle for lowest_l_trajectories.dat.

static double **g_lowestL_L_buf = NULL

Buffer for angular momentum of low-L particles [particle_idx][timestep_in_buffer].

static double **g_lowestL_r_buf = NULL

Buffer for radius of low-L particles [particle_idx][timestep_in_buffer].

static unsigned long int g_master_seed = 0

Master seed for the simulation, if provided.

static int g_master_seed_provided = 0

Flag: 1 if the --master-seed option is provided.

static double **g_mu_buf = NULL

Buffer for radial direction cosine mu [particle_idx][timestep_in_buffer].

static double g_nfw_profile_falloff_factor = FALLOFF_FACTOR_NFW_DEFAULT

Falloff concentration parameter \(C\) for the NFW profile. Populated from g_falloff_factor_param.

static double g_nfw_profile_halo_mass = HALO_MASS_NFW

Total halo mass for the NFW profile ( \(M_{\odot}\)). Populated from g_halo_mass_param.

static double g_nfw_profile_rc = RC_NFW_DEFAULT

Scale radius \(r_c\) for the NFW profile (kpc). Populated from g_scale_radius_param.

static double g_nfw_profile_rmax_norm_factor = CUTOFF_FACTOR_NFW_DEFAULT

Cutoff radius factor for the NFW profile. Populated from g_cutoff_factor_param.

static int g_om_aniso_factor_provided = 0

Flag: 1 if the --aniso-factor or --aniso-betascale option is provided.

static double g_om_anisotropy_radius = 0.0

Physical anisotropy radius \(r_a\) (kpc)

static double g_om_ra_scale_factor = 0.0

Anisotropy radius as multiple of scale radius.

static int *g_particle_scatter_state = NULL

Tracks recent scatter history for Adams-Bashforth integrator state reset. Indexed by original particle ID. 0=normal AB, 1=scattered in the previous step (use AB1-like step), 2=one step has passed since scattering (use AB2-like step).

static double g_precalc_force_const = 0.0

Pre-calculated gravitational force constant.

Combines -G * M_total / npts * unit conversions to avoid repeated division in force calculations. Updated once before main time loop.

static char g_profile_type_str[16] = "nfw"

Profile type string ("nfw", "cored", or "hernquist"), default "nfw".

static int g_profile_type_str_provided = 0

Flag: 1 if the --profile option is provided.

static int g_quadsort_wins = 0

Number of times quadsort was faster.

static int g_radix_wins = 0

Number of times radix sort was faster.

char *g_restart_file_override = NULL

Optional explicit restart file path (--restart-file).

int g_restart_initial_timestep = 0

Actual timestep number from which restart continues (e.g., 10001).

int g_restart_mode_active = 0

Flag: 1 when actively restarting simulation, 0 when only checking or not restarting.

int g_restart_snapshots_is_count = 0

Flag: 1 if restart_completed_snapshots represents a count of snapshots; 0 if it is an index.

static gsl_rng *g_rng = NULL

GSL Random Number Generator state for IC generation.

static double g_scale_radius_param = RC

Generalized scale radius (kpc), defaults to RC macro.

static int g_scale_radius_param_provided = 0

Flag: 1 if the --scale-radius option is provided.

int g_sidm_execution_mode = 1

SIDM execution mode: 0 for serial, 1 for parallel (default).

static double g_sidm_kappa = 50.0

SIDM opacity kappa (cm \(^2\)/g), default 50.0.

static int g_sidm_kappa_provided = 0

Flag: 1 if the --sidm-kappa option is provided.

int g_sidm_max_interaction_range = 10

Maximum number of neighbors to check for SIDM scattering. Default is 10.

static unsigned long int g_sidm_seed = 0

Seed used for SIDM calculations.

static const char *g_sidm_seed_filename_base = "data/last_sidm_seed"

Base name for SIDM seed file.

static int g_sidm_seed_provided = 0

Flag: 1 if the --sidm-seed option is provided.

static FILE *g_single_traj_file = NULL

File handle for single_trajectory.dat.

static int g_sort_call_count = 0

Total number of sort calls made.

static double **g_sort_columns_buffer = NULL

Global persistent buffer for particle data transposition during sorting.

static int g_sort_columns_buffer_npts = 0

Stores the number of particles for which the persistent buffer is allocated.

The buffer is reallocated if the number of simulation particles changes.

static double *g_time_buf = NULL

Buffer for time values [timestep_in_buffer].

static double *g_time_snapshots = NULL

Arrays for tracking total system energy diagnostics over time.

Simulation time at each diagnostic snapshot (Myr).

static double *g_total_E = NULL

Total energy (KE + PE) of the system at each snapshot.

static double g_total_insertion_time = 0.0

Total time spent in insertion benchmarks (ms)

static double *g_total_KE = NULL

Total kinetic energy of the system at each snapshot.

static double *g_total_PE = NULL

Total potential energy of the system at each snapshot.

static double g_total_quadsort_time = 0.0

Total time spent in quadsort benchmarks (ms)

static double g_total_radix_time = 0.0

Total time spent in radix benchmarks (ms)

long long g_total_sidm_scatters = 0

Global counter for total SIDM scatters.

static FILE *g_traj_file = NULL

File handle for trajectories.dat, kept open for performance.

static double **g_trajectories_buf = NULL

Buffer for radius [particle_idx][timestep_in_buffer].

static int g_trajectory_buffer_index = 0

Current index for writing into the trajectory buffers.

static int g_trajectory_buffer_size = 0

Size of the trajectory buffers (in timesteps), synchronized with snapshot buffer size.

int g_use_graph_coloring_sidm = 0

Use graph coloring algorithm for parallel SIDM (eliminates double-booking).

static int g_use_hernquist_aniso_profile = 0

Flag to use anisotropic Hernquist profile for ICs.

static int g_use_hernquist_numerical = 0

Flag to use numerical Hernquist profile (OM-compatible).

static int g_use_nfw_profile = 0

Flag to select NFW profile for ICs (1=NFW, 0=use other profile flags).

static int g_use_numerical_isotropic = 0

Flag for numerical Hernquist with true f(E).

static int g_use_om_profile = 0

Osipkov-Merritt Anisotropy Model Implementation.

The Osipkov-Merritt (OM) model provides radially-varying velocity anisotropy characterized by the anisotropy parameter:

\(\beta(r) = r^2/(r^2 + r_a^2)\)

where \(r_a\) is the anisotropy radius. The model transitions from isotropic ( \(\beta=0\)) at \(r=0\) to radially-biased ( \(\beta \rightarrow 1\)) as \(r \rightarrow \infty\).

Implementation uses the augmented density method:

• Define augmented density: \(\rho_Q(r) = \rho(r)(1 + r^2/r_a^2)\)

• Use \(Q = E - L^2/(2r_a^2)\) as the OM invariant

• Apply Eddington inversion to \(\rho_Q\) to get \(f(Q)\)

• Sample velocities in pseudo-space where Q-surfaces are spheres

• Transform to physical velocities using mapping: \(v_r = w \cos(\theta)\), \(v_t = w \sin(\theta)/\sqrt{1 + r^2/r_a^2}\)

Note

Compatible with NFW, Cored, and numerical Hernquist profiles. Flag: 1 if OM model is requested

static double **g_velocities_buf = NULL

Buffer for radial velocity [particle_idx][timestep_in_buffer].

static int *ID_block = NULL

Particle ID block (original IDs).

static float *L_block = NULL

Global arrays for particle data processing (block storage).

Angular momentum block.

static double Lcompare = 0.05

Reference \(L\) value for closest-match mode (Mode 1).

static int nlowest = 5

Variables for tracking low angular momentum particles.

Number of lowest angular momentum particles to track.

static double normalization

Mass normalization factor for energy calculations. Calculated based on the integral of the density profile.

static int num_traj_particles = 10

Number of particles to track in trajectories.dat (by original ID).

static const int PARALLEL_SORT_DEFAULT_SECTIONS = 8

Default number of sections when OpenMP is unavailable or reports few threads.

Provides a baseline level of partitioning even in limited thread environments.

static const int PARALLEL_SORT_MIN_CHUNK_SIZE_THRESHOLD = 128

Minimum average chunk size threshold for parallel sort operation.

If the calculated size of each chunk (total elements / number of sections) falls below this threshold, the algorithm reverts to serial sorting to avoid overhead from managing very small parallel tasks.

static const int PARALLEL_SORT_MIN_CORRECTNESS_OVERLAP = 32

Minimum required overlap between adjacent sort sections (in elements).

This ensures sufficient overlap for correct merging of sorted sections, particularly for sparse data distributions or when the calculated proportional overlap is small.

static const int PARALLEL_SORT_OVERLAP_DIVISOR = 8

Divisor for calculating proportional overlap between sort sections.

Overlap is calculated as chunk_size / OVERLAP_DIVISOR to scale with data size.

static const int PARALLEL_SORT_SECTIONS_PER_THREAD = 2

Number of sort sections per OpenMP thread for workload distribution.

Multiplier used to determine total section count from available threads.

static float *phi_block = NULL

Phi angle block.

static float *R_block = NULL

Radial position block.

static int *Rank_block = NULL

Particle rank (sorted position) block.

static const char *readInitFilename = NULL

Filename to read initial conditions from (if doReadInit=1).

static int *scatter_count_block = NULL

Scatter count block.

int skip_file_writes = 0

Skip file writes during simulation restart.

static int use_closest_to_Lcompare = 0

Mode selector (0=lowest \(L\), 1=closest to \(L_{compare}\)). Default is lowest \(L\).

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 float *Vrad_block = NULL

Radial velocity block.

static const char *writeInitFilename = NULL

Filename to write initial conditions to (if doWriteInit=1).

struct cored_potential_params

Public Members

double M

double rc

struct fE_integrand_params_hernquist_t

Parameters for the Hernquist distribution function integrand fEintegrand_hernquist.

This structure passes all necessary data, including pre-computed splines for \(r(\Psi)\) and \(M(r)\), physical constants, and profile-specific parameters, to the GSL integration routine for calculating \(I(E)\) or \(I(Q)\).

Public Members

gsl_interp_accel *accel_M_of_r

Accelerator for the M(r) spline.

gsl_interp_accel *accel_r_of_Psi

Accelerator for the \(r(-\Psi)\) spline.

double const_G_universal

Universal gravitational constant G.

double E_current_shell

Energy E (isotropic) or Q (OM) of the current shell for which I is being computed.

double hernquist_a_scale

Scale radius \(a\) of the Hernquist profile.

double hernquist_normalization

Density normalization constant for the Hernquist profile.

double Psimax_global

Maximum potential value for physical range validation.

double Psimin_global

Minimum potential value for physical range validation.

gsl_spline *spline_M_of_r

Spline for M(r), enclosed mass as a function of radius.

gsl_spline *spline_r_of_Psi

Spline for \(r(-\Psi)\): radius as function of negated potential (negation ensures monotonicity for GSL).

int use_om

Flag: 1 if using OM (E_current_shell represents Q), 0 if isotropic.

struct fE_integrand_params_NFW_t

Parameters for the NFW distribution function integrand fEintegrand_nfw.

This structure passes all necessary data, including pre-computed splines for \(r(\Psi)\) and \(M(r)\), physical constants, and profile-specific parameters, to the GSL integration routine for calculating \(I(E)\).

Public Members

gsl_interp_accel *accel_M_of_r

Accelerator for the M(r) spline.

gsl_interp_accel *accel_r_of_Psi

Accelerator for the \(r(-\Psi)\) spline.

double const_G_universal

Universal gravitational constant G.

double E_current_shell

Energy E (isotropic) or Q (OM) of the current shell for which I is being computed.

double profile_falloff_C_const

Falloff concentration parameter \(C\) for power-law cutoff in NFW profile.

double profile_nt_norm_const

Density normalization constant (nt_nfw) for the NFW profile.

double profile_rc_const

Scale radius \(r_c\) of the NFW profile.

double Psimax_global

Maximum potential value for physical range validation.

double Psimin_global

Minimum potential value for physical range validation.

gsl_spline *spline_M_of_r

Spline for M(r), enclosed mass as a function of radius.

gsl_spline *spline_r_of_Psi

Spline for \(r(-\Psi)\): radius as function of negated potential (negation ensures monotonicity for GSL).

int use_om

Flag: 1 if using OM (E_current_shell represents Q), 0 if isotropic.

struct fEintegrand_params

Parameters for energy integration calculations.

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

Public Members

double E

Energy value \(E\) for isotropic or \(Q\) value for OM.

gsl_interp_accel *massarray

Accelerator for mass lookups \(M(r)\).

gsl_interp_accel *rofPsiarray

Accelerator for radius lookup from potential \(r(\Psi)\).

gsl_spline *splinemass

Interpolation spline for enclosed mass \(M(r)\).

gsl_spline *splinePsi

Interpolation spline for potential \(\Psi(r)\).

int use_om

Flag: 1 if using OM (E represents Q), 0 if isotropic (E represents energy).

struct hernquist_params

Struct to pass parameters to the GSL Hernquist density integrand.

Public Members

double a

Scale Radius.

double M

Total Mass.

struct hernquist_potential_params

Public Members

double a

double M

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

int idx

Original particle index (before sorting by L).

double L

Angular momentum value (or squared difference from Lcompare).

int sign

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

struct nfw_potential_params

Public Members

double falloff

double M

double normalization

double Phi0

double rs

struct om_mu_integral_params

Parameters for OM \(\mu\)-integral in df_fixed_radius output.

Public Members

gsl_interp_accel *fofQ_accel

gsl_interp *fofQ_interp

double *I_array

double Psi_rf_codes

double Psimax

double Psimin

double *Q_array

double r_over_ra_sq

double v_sq

struct PartData

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

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

float angmom

int original_index

float rad

int rank

float vrad

struct Psiintegrand_params

Parameters for Psiintegrand to support profile-specific mass integrands.

Allows Psiintegrand to call the appropriate mass integrand function based on the selected density profile.

Public Members

double (*massintegrand_func)(double, void*)

Function pointer to profile-specific mass integrand.

void *params_for_massintegrand

Parameters for the mass integrand function.

struct RrPsiPair

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

Orchestrates the overall simulation workflow:

1. Parses command-line arguments using a flexible option system

2. Sets up initial conditions (particle distribution, system parameters)

3. Executes the selected integration method for trajectory evolution

4. Outputs data products (particle states, phase diagrams, energy tracking)

5. Handles restart/resume functionality when requested

Note

This is a highly parallelized application that takes advantage of OpenMP when available. Performance scales with the number of available cores.

Warning

Some simulation configurations can be very memory-intensive. For large particle counts ( \(>10^6\)), ensure sufficient RAM is available.

Param argc:

[in] Standard argument count from the command line.

Param argv:

[in] Standard array of argument strings from the command line.

Return:

Exit code: 0 for successful execution, non-zero for errors.

Public Members

double psi

Corresponding potential value or y-axis value.

double rr

Radius value or x-axis value for sorting.

struct threevector

Three-dimensional vector structure for particle physics calculations.

Represents velocity and position vectors in 3D space for self-interacting dark matter simulations. Components are stored in Cartesian coordinates.

Public Members

double x

x-component of the vector.

double y

y-component of the vector.

double z

z-component of the vector.