Metadata-Version: 2.4
Name: mcfast
Version: 0.1.5
Classifier: Programming Language :: Rust
Classifier: Programming Language :: Python :: Implementation :: CPython
Classifier: Programming Language :: Python :: Implementation :: PyPy
License-File: LICENSE.md
Summary: A set of utilities and Rust-accelerated functions for the McFACTS AGN disk simulation project
Keywords: astronomy,physics
Author-email: Nicolas Roman Posner <nicolasposner@gmail.com>
Maintainer-email: Nicolas Roman Posner <nicolasposner@gmail.com>
License: MIT
Requires-Python: >=3.8
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM

# McFAST v0.1.5

Utilities and accelerated functions designed for use by the McFACTS team at CUNY.

Still under early alpha development. Modules and functions that are not explicitly documented here are not yet ready for use.

## Tau Incline/Eccentricity

In order to accelerate the tau_inc_dyn and tau_ecc_dyn functions in `mcfacts/physics/disk_capture.py`, we provide `tau_inc_dyn_helper` and `tau_ecc_dyn_helper` functions that perform the vectorized array calculations in an optimized Rust function. The example implementation for `tau_inc_dyn_optimized` follows: it has the same call signature, takes the same inputs, and provides the same output as `tau_inc_dyn`, but is 2-3x faster.

This function and its `ecc` variant can be further optimized by speeding up or eliminating the calls to si_from_r_g, which now make up the bulk of the runtime in each function (approximately 2/3rds).

In particular, the `tau_semi_lat` functionality in `tau_ecc_dyn` is sped up at least 120x, with tau_semi_lat initially being timed at 6 seconds, while the `tau_ecc_dyn_helper` functionality that subsumes it is timed at 0.05 seconds.

```python
def tau_inc_dyn_optimized(smbh_mass, disk_bh_retro_orbs_a, disk_bh_retro_masses, omega, ecc, inc, disk_surf_density_func, r_g_in_meters):
    """Computes inclination damping timescale from actual variables; used only for scaling.
    Uses Rust-accelerated helper functions for the calculation: compare to tau_inc_dyn
    """
    # throw most things into SI units (that's right, ENGINEER UNITS!)
    #    or more locally convenient variable names
    SI_smbh_mass = smbh_mass * u.Msun.to("kg")  # kg
    SI_semi_maj_axis = si_from_r_g(smbh_mass, disk_bh_retro_orbs_a, r_g_defined=r_g_in_meters).to("m").value
    SI_orbiter_mass = disk_bh_retro_masses * u.Msun.to("kg")  # kg
    cos_omega = np.cos(omega)

    disk_surf_res = disk_surf_density_func(disk_bh_retro_orbs_a)

    tau_i_dyn = tau_inc_dyn_helper(SI_smbh_mass, SI_orbiter_mass, ecc, inc, cos_omega, disk_surf_res, SI_semi_maj_axis)

    assert np.isfinite(tau_i_dyn).all(), \
        "Finite check failure: tau_i_dyn"

    return tau_i_dyn
```

## Merge Tree
For M. McCarthy's black hole merge tracking feature.

```python
import mcfast
from mcfast import MergeForest

# creates a forest of all black holes
forest = MergeForest("./data/", "galaxy_state_*")

bh_uuid: str = "2b422064-0a84-4687-a542-395dcb61cd4f"

# get the immediate child of a given UUID, if the UUID is valid
descendant_uuid = forest.get_descendant(bh_uuid)

# get the immediate parents of a given UUID, if the UUID is valid
(parent1_uuid, parent2_uuid) = forest.get_parents(bh_uuid)

# get the full ancestry list of a given UUID (in DFS order), 
# assuming the UUID has ancestors
ancestor_list = forest.get_ancestors(bh_uuid)


# get the full list of black holes between the UUID and 
# the final product (root)
descent_path = forest.get_lineage_to_root(bh_uuid)

# get the generation of a given UUID, where an initialized BH is generation 0
# note: implementation doesn't fully line up with existing
# reference implementation, under review
generation = forest.get_generation(bh_uuid)

# get all 'root' black holes (products remaining at the final tick)
roots = forest.roots()

# get all 'leaf' black holes (initialized BHs with no ancestors)
leaves = forest.leaves()

# get all 'singleton' black holes (neither ancestors nor children)
singletons = forest.singletons()

if forest.contains(bh_uuid):
    print("It exists!")
else:
    print("It doesn't exist!")

# get total number of nodes in the 'forest'
total = len(forest)

# serialize and deserialize with pickle
import pickle
 
# Save now
with open("forest.pkl", "wb") as f:
    pickle.dump(forest, f)

# Load later
with open("forest.pkl", "rb") as f:
    restored = pickle.load(f)

# sample: 330-node forest serializes down to 16kb
```


