Quick Start

A very simple example in both Python and C++. The galaxy_interaction.hdf5 snapshot comes from the public Arepo test suite. For more detailed examples, see the Tutorials.

Load data and project

Python

import numpy as np
import h5py as h5
import vortrace as vt

# Load an Arepo snapshot
f = h5.File("snapshot.hdf5", "r")
pos = np.array(f["PartType0/Coordinates"])
rho = np.array(f["PartType0/Density"])
mass = np.array(f["PartType0/Masses"])
vol = mass / rho
BoxSize = f["Parameters"].attrs["BoxSize"]
f.close()

# Create a projection cloud
pc = vt.ProjectionCloud(
    pos, rho, vol=vol,
    boundbox=[0, BoxSize, 0, BoxSize, 0, BoxSize],
)

# Make a grid projection
L = 75.0
extent = [BoxSize / 2 - L / 2, BoxSize / 2 + L / 2]
bounds = [0, BoxSize]
npix = 256

image = pc.grid_projection(extent, npix, bounds, center=None)

# Plot the result
fig, ax, im = vt.plot.plot_grid(
    image,
    extent=[-L / 2, L / 2, -L / 2, L / 2],
    label="Column density",
)

C++

#include <vortrace/vortrace.hpp>
#include <vector>
#include <cmath>

int main() {
    // --- Load your particle data into flat arrays ---
    // pos:    [x0, y0, z0, x1, y1, z1, ...]
    // fields: [rho0, rho1, ...]
    size_t npart = /* number of particles */;
    std::vector<double> pos(3 * npart);
    std::vector<double> fields(npart);
    // ... fill pos and fields from your data source ...

    double BoxSize = /* your box size */;

    // Build a point cloud and kD-tree
    std::array<double, 6> subbox = {0, BoxSize, 0, BoxSize, 0, BoxSize};
    PointCloud cloud;
    cloud.loadPoints(pos.data(), npart, fields.data(), npart, 1, subbox);
    cloud.buildTree();

    // Set up a grid of rays
    double L = 75.0;
    double cx = BoxSize / 2, cy = BoxSize / 2;
    size_t npix = 256;
    double dx = L / npix;

    std::vector<Float> starts(3 * npix * npix);
    std::vector<Float> ends(3 * npix * npix);

    for (size_t iy = 0; iy < npix; iy++) {
        for (size_t ix = 0; ix < npix; ix++) {
            size_t idx = 3 * (iy * npix + ix);
            double x = cx - L / 2 + (ix + 0.5) * dx;
            double y = cy - L / 2 + (iy + 0.5) * dx;
            starts[idx]     = x;
            starts[idx + 1] = y;
            starts[idx + 2] = 0.0;
            ends[idx]       = x;
            ends[idx + 1]   = y;
            ends[idx + 2]   = BoxSize;
        }
    }

    // Run the projection
    Projection proj(starts.data(), ends.data(), npix * npix);
    proj.makeProjection(cloud, ReductionMode::Sum);

    // Access the result
    const auto& data = proj.getProjectionData();
    // data[i] is the integrated column density for ray i
}

Result