# matSpringDampMass

Builds the system matrices for the granular eigenvalue problem
`M*x''(t) + Gamma*x'(t) + K*x(t) = 0`: the stiffness matrix (Hessian) `matSpring`,
the damping matrix `matDamp`, and the mass matrix `matMass` for a 2D or 3D
packing with Hookean contacts. Degrees of freedom are ordered
`[x1, y1, (z1), x2, y2, (z2), ...]`, so the matrices are `(dof x dof)` with
`dof = N*dim`.

## Default Syntax

```matlab
packing = data.packing;     % struct, loaded from a pack.m results .mat file
[K, Gamma, M] = matSpringDampMass(packing, dampingConstant, springConstant, mass)
```

With `dampingConstant = 1`, `springConstant = 1`, `mass = 1` (all defaults).
The packing struct is detected automatically: 2D vs 3D from the presence of
`vecPosZ`, new format (`vecPosX`/`scalBoxWidthX`/...) vs old format via the
`oldPackingFormat` option.

## Options

```matlab
opts.algorithm        = "cell";     % "cell" (default, linear scaling, fast for large N)
                                    % "particle" (all pairs, quadratic scaling, fine for small N)
opts.periodic         = false;      % logical, periodic boundary conditions
opts.sparseOutput     = false;      % logical, return sparse matrices
opts.uniformMass      = false;      % logical, all particles get the same mass
opts.oldPackingFormat = false;      % logical, packing uses x/y(/z), Dn, Lx/Ly(/Lz) fields
[K, Gamma, M] = matSpringDampMass(packing, dampingConstant, springConstant, mass, opts)
```

- `algorithm = "cell"` buckets particles into cells of width `4*max(radii)` and
  only checks the neighboring cells — use for large packings. `"particle"`
  checks every pair — use for small ones.
- `periodic` must be true for the cell algorithm to run (the non-periodic cell
  branch is not implemented). For non-periodic 2D packings, wall contacts
  (particle within its radius of a wall) add `K` to the diagonal.
- `sparseOutput` returns `sparse` K/Gamma/M — needed before passing to
  `polyeig` for large systems.
- `uniformMass` gives every particle mass `mass`; otherwise mass =
  `mass * (unit-ball volume) * r^dim`.

## Examples

```matlab
data = load('data/eigenData/results_2D_iso_N100_P0.1_Seed5_gamma_3.00e-03.mat');

% Damped system matrices, dense
[K, Gamma, M] = matSpringDampMass(data.packing, 0.01, 100, 1);

% Sparse, periodic — ready for polyeig
opts.periodic      = true;
opts.sparseOutput  = true;
[K, Gamma, M] = matSpringDampMass(data.packing, 0.001, 100, 1, opts);

[eigenVectors, eigenValues] = polyeig(K, Gamma, M);

% Old-format packing struct (x/y/Dn/Lx/Ly)
packing = struct('x', x, 'y', y, 'Dn', Dn, 'Lx', Lx, 'Ly', Ly);
opts.oldPackingFormat = true;
opts.periodic         = true;
[K, Gamma, M] = matSpringDampMass(packing, 1, 100, 1, opts);

% Small packing, all-pairs check
opts2 = struct('algorithm', "particle");
[K, Gamma, M] = matSpringDampMass(packing, 1, 100, 1, opts2);
```