diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index b60baa7..4d798e6 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -5,27 +5,24 @@ on: [pull_request, push]
 jobs:
   test:
     runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        python: ['3.10', '3.11', '3.12']
     steps:
       - uses: actions/checkout@v4
 
-      - uses: actions/setup-python@v5
+      - uses: conda-incubator/setup-miniconda@v3
         with:
-          python-version: ${{ matrix.python }}
-
-      - uses: julia-actions/setup-julia@v2
-        with:
-          version: '1.10'
-
-      - uses: julia-actions/cache@v2
+          environment-file: environment.yml
+          auto-activate-base: false
+          miniforge-version: latest
+          use-mamba: true
 
       - name: Install MAGEMin_C
-        run: julia -e 'using Pkg; Pkg.add("MAGEMin_C")'
+        shell: bash -l {0}
+        run: phasetools-julia-setup --install
 
-      - name: Install phasetools
-        run: python -m pip install -e .
+      - name: Verify environment
+        shell: bash -l {0}
+        run: phasetools-julia-setup --check
 
       - name: Run tests
+        shell: bash -l {0}
         run: python -m unittest discover tests -v
diff --git a/AGENTS.md b/AGENTS.md
index ef94c3a..c509106 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -2,10 +2,19 @@
 
 ## Quick start
 
+**Conda / mamba (Linux & CI):**
+```bash
+mamba env create -f environment.yml
+conda activate phasetools
+phasetools-julia-setup --install
+phasetools-julia-setup --check
+```
+
+**pip + juliaup (cross-platform, incl. macOS):**
 ```bash
 python -m pip install -e .
-phasetools-julia-setup --check    # verify Julia + MAGEMin_C
 phasetools-julia-setup --install  # install MAGEMin_C in Julia
+phasetools-julia-setup --check    # verify Julia + MAGEMin_C
 python3 -m unittest discover tests   # run all tests (mock-based, no Julia needed)
 ```
 
@@ -282,9 +291,10 @@ The `ox.count("O")` bug undercounts oxygen, inflating the6-O scale factor by ~60
 - Files: `test_redox_logic.py`, `test_site_occupancy.py`, `test_fe_basis.py`, `test_fractionation.py`, `test_iron_oxide_conversions.py`, `test_solvus_instances.py`, `test_lmo_fix.py`.
 - New public functions must be documented in the **directory's README.md** (e.g., `calculators/README.md`, `core/README.md`) **and** the package structure above must be updated.
 
-## Requirements (from `setup.py`)
+## Requirements (from `pyproject.toml`)
 
 `pandas`, `numpy`, `matplotlib`, `scipy`, `molmass`, `juliacall`. Python ≥ 3.10.
+Extras: `[jupyter]` (jupyterlab, ipykernel, nbconvert), `[dev]` (ruff).
 
 ## Tutorials
 
diff --git a/README.md b/README.md
index f747d1a..4c6dbf9 100644
--- a/README.md
+++ b/README.md
@@ -25,24 +25,29 @@
 
 ## 📦 Installation
 
-### 1. Python Package
-Install `phasetools` in editable mode from the repository root:
+### 1. Conda / mamba (Linux & CI)
 
 ```bash
-python -m pip install -e .
+mamba env create -f environment.yml
+conda activate phasetools
+phasetools-julia-setup --install
+phasetools-julia-setup --check
 ```
 
-### 2. Julia Environment Setup
-`phasetools` requires Julia and the `MAGEMin_C` package. Use the built-in helper to check your environment:
+*Note: On macOS (especially Apple Silicon), conda-forge's `julia` package is unavailable — use the pip + juliaup path below instead.*
+
+### 2. pip + juliaup (cross-platform, incl. macOS)
 
 ```bash
+python -m pip install -e .
 phasetools-julia-setup --check
+phasetools-julia-setup --install
 ```
 
-If `MAGEMin_C` is missing, install it automatically:
+To also install Jupyter (for tutorials) and dev tooling (ruff):
 
 ```bash
-phasetools-julia-setup --install
+python -m pip install -e ".[jupyter,dev]"
 ```
 
 *Note: If Julia is not installed on your system, please download it from [julialang.org](https://julialang.org/downloads/).*
diff --git a/Tutorials/Garnet/1-Functionality_overview.ipynb b/Tutorials/Garnet/1-Functionality_overview.ipynb
index cf2f9f1..11a1651 100644
--- a/Tutorials/Garnet/1-Functionality_overview.ipynb
+++ b/Tutorials/Garnet/1-Functionality_overview.ipynb
@@ -26,9 +26,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from phasetools import MAGEMinGarnetCalculator\n",
-    "import juliacall\n",
-    "import matplotlib.pyplot as plt"
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from phasetools import MAGEMinGarnetCalculator"
    ]
   },
   {
@@ -207,6 +207,7 @@
    "outputs": [],
    "source": [
     "from phasetools.utils import create_PTt_path\n",
+    "\n",
     "P = [5., 6., 7., 8., 9., 10.]  # kbar\n",
     "T = [500., 600., 700., 800., 900., 1000.]  # °C\n",
     "t = [1., 2., 3., 4., 5., 6.]  # Myr\n",
@@ -290,7 +291,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from juliacall import Main as jl, convert as jlconvert\n",
+    "from juliacall import Main as jl\n",
+    "from juliacall import convert as jlconvert\n",
+    "\n",
     "rm_j = jlconvert(jl.Vector[jl.String], ['chl'])\n",
     "rmp_MM = MAGEMin_C.remove_phases(rm_j, 'mp')"
    ]
diff --git a/Tutorials/Garnet/2-Overlapping_isopleths.ipynb b/Tutorials/Garnet/2-Overlapping_isopleths.ipynb
index 0805319..e6eefdc 100644
--- a/Tutorials/Garnet/2-Overlapping_isopleths.ipynb
+++ b/Tutorials/Garnet/2-Overlapping_isopleths.ipynb
@@ -28,15 +28,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import pandas as pd\n",
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt\n",
-    "from phasetools import PhasePTEstimator\n",
-    "from phasetools.utils import *\n",
-    "import phasetools\n",
     "import os\n",
     "\n",
-    "from juliacall import Main as jl, convert as jlconvert"
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "from phasetools import PhasePTEstimator\n",
+    "from phasetools.utils import *"
    ]
   },
   {
@@ -122,6 +121,7 @@
    "outputs": [],
    "source": [
     "from phasetools import MAGEMinGarnetCalculator\n",
+    "\n",
     "garnet_calculator = MAGEMinGarnetCalculator(db=db, dataset=636)\n",
     "garnet_calculator.setup_bulk_composition(Xoxides, X, sys_in=sys_in)\n",
     "\n",
diff --git a/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb b/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb
index 2894c8a..08a364d 100644
--- a/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb
+++ b/Tutorials/Garnet/3-Garnet_cpx_thermometry.ipynb
@@ -52,13 +52,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
     "\n",
-    "from phasetools import MAGEMinPTGridCalculator\n",
-    "\n",
-    "from phasetools import MAGEMin_C\n",
-    "from juliacall import Main as jl, convert as jlconvert"
+    "from phasetools import MAGEMinPTGridCalculator"
    ]
   },
   {
@@ -569,7 +566,7 @@
    "source": [
     "from phasetools.core.phase_properties import get_phase_fe_split\n",
     "\n",
-    "for i in range(0, 5):\n",
+    "for i in range(5):\n",
     "      i0 = int(np.flatnonzero(valid)[i])  # first valid grid point\n",
     "      o0 = out[i0]\n",
     "      cpx_split = get_phase_fe_split(o0, cpx)\n",
diff --git a/Tutorials/General/1-Getting_started.ipynb b/Tutorials/General/1-Getting_started.ipynb
index 43d5782..733a57c 100644
--- a/Tutorials/General/1-Getting_started.ipynb
+++ b/Tutorials/General/1-Getting_started.ipynb
@@ -24,10 +24,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
-    "import juliacall\n",
-    "from juliacall import Main as jl, convert as jlconvert\n",
-    "import phasetools\n",
+    "from juliacall import Main as jl\n",
+    "from juliacall import convert as jlconvert\n",
+    "\n",
     "from phasetools import MAGEMin_C"
    ]
   },
diff --git a/Tutorials/General/2-Determining_fractions.ipynb b/Tutorials/General/2-Determining_fractions.ipynb
index 4f871af..fce609a 100644
--- a/Tutorials/General/2-Determining_fractions.ipynb
+++ b/Tutorials/General/2-Determining_fractions.ipynb
@@ -28,8 +28,7 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
-    "import juliacall\n",
-    "from juliacall import Main as jl, convert as jlconvert\n",
+    "\n",
     "import phasetools\n",
     "from phasetools import MAGEMin_C"
    ]
@@ -103,6 +102,7 @@
    "outputs": [],
    "source": [
     "import matplotlib.pyplot as plt\n",
+    "\n",
     "plt.plot(temp, liq_frac_vals, c='k', ls=':')\n",
     "plt.scatter(liquidus_T, 1, marker='x', c='red', label='Liquidus')\n",
     "plt.scatter(solidus_T, 0, marker='x', c='green', label='Solidus')\n",
@@ -192,6 +192,7 @@
    "source": [
     "from scipy.optimize import root_scalar\n",
     "\n",
+    "\n",
     "def melt_fraction(T):\n",
     "    out = MAGEMin_C.single_point_minimization(P, T, data)\n",
     "    return phasetools.phase_frac(phase=\"liq\", MAGEMinOutput=out, sys_in='mol')\n",
@@ -241,6 +242,7 @@
    "outputs": [],
    "source": [
     "from phasetools import PhaseFunctions\n",
+    "\n",
     "phaseFunctions = PhaseFunctions(db=db, dataset=636)\n",
     "phaseFunctions.setup_bulk_composition(list(out.oxides), list(out.bulk), sys_in='mol')"
    ]
diff --git a/Tutorials/General/3-Lunar_Magma_Ocean.ipynb b/Tutorials/General/3-Lunar_Magma_Ocean.ipynb
index 5167c25..ec07f5f 100644
--- a/Tutorials/General/3-Lunar_Magma_Ocean.ipynb
+++ b/Tutorials/General/3-Lunar_Magma_Ocean.ipynb
@@ -24,10 +24,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
     "import matplotlib.pyplot as plt\n",
-    "from phasetools.models.magma_ocean import MagmaOcean\n",
-    "import sys"
+    "\n",
+    "from phasetools.models.magma_ocean import MagmaOcean"
    ]
   },
   {
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 0000000..62ce049
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,10 @@
+name: phasetools
+channels:
+  - conda-forge
+  - nodefaults
+dependencies:
+  - python=3.11
+  - julia
+  - pip
+  - pip:
+      - -e ".[jupyter,dev]"
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..bb7c118
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,45 @@
+[build-system]
+requires = ["setuptools>=61.0", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "phasetools"
+version = "0.1.0"
+description = "A python package to perform MAGEMin calculations in python"
+readme = "README.md"
+license = { file = "LICENSE" }
+requires-python = ">=3.10"
+authors = [{ name = "Ben Knight", email = "ben.knight@curtin.edu.au" }]
+classifiers = [
+    "Development Status :: 2 - Pre-Alpha",
+    "Intended Audience :: Developers",
+    "License :: OSI Approved :: GNU Lesser General Public License v3 or later (LGPLv3+)",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.10",
+    "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
+    "Topic :: Software Development :: Libraries :: Python Modules",
+    "Topic :: Scientific/Engineering",
+    "Natural Language :: English",
+]
+dependencies = [
+    "pandas>=1.0.0",
+    "numpy>=1.18.0",
+    "matplotlib>=3.0.0",
+    "scipy>=1.4.0",
+    "molmass>=2024.0.0",
+    "juliacall>=0.9",
+]
+
+[project.optional-dependencies]
+jupyter = ["jupyterlab", "ipykernel", "nbconvert"]
+dev = ["ruff"]
+
+[project.scripts]
+phasetools-julia-setup = "phasetools.julia_setup:main"
+
+[tool.setuptools.packages.find]
+where = ["src"]
+
+[tool.ruff]
+extend-exclude = ["Tutorials"]
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 8c0c38c..0000000
--- a/setup.py
+++ /dev/null
@@ -1,38 +0,0 @@
-from setuptools import setup, find_packages
-
-setup(
-    name='phasetools',
-    version='0.1.0',
-    author='Ben Knight',
-    author_email='ben.knight@curtin.edu.au',
-    description='A python package to perform MAGEMin calculations in python',
-    packages=find_packages(where='src'),
-    package_dir={'': 'src'},
-    install_requires=[
-        'pandas>=1.0.0',
-        'numpy>=1.18.0',
-        'matplotlib>=3.0.0',
-        'scipy>=1.4.0',
-        'molmass>=2024.0.0',
-        'juliacall>=0.9',
-    ],
-    entry_points={
-        'console_scripts': [
-            'phasetools-julia-setup=phasetools.julia_setup:main',
-        ],
-    },
-    license='LGPL-3.0-or-later',
-    classifiers=[
-        'Development Status :: 2 - Pre-Alpha',
-        'Intended Audience :: Developers',
-        'License :: OSI Approved :: GNU Lesser General Public License v3 or later (LGPLv3+)',
-        'Programming Language :: Python :: 3',
-        'Programming Language :: Python :: 3.10',
-        'Programming Language :: Python :: 3.11',
-        'Programming Language :: Python :: 3.12',
-        'Topic :: Software Development :: Libraries :: Python Modules',
-        'Topic :: Scientific/Engineering',
-        'Natural Language :: English',
-    ],
-    python_requires='>=3.10',
-)
\ No newline at end of file
diff --git a/src/phasetools/__init__.py b/src/phasetools/__init__.py
index 52357ac..0ebb854 100644
--- a/src/phasetools/__init__.py
+++ b/src/phasetools/__init__.py
@@ -3,11 +3,12 @@
 ### Import juliacall to access MAGEMin_C
 try:
     import juliacall
-    from juliacall import Main as jl, convert as jlconvert
+    from juliacall import Main as jl
+    from juliacall import convert as jlconvert
 
     MAGEMin_C = juliacall.newmodule("MAGEMin_C")
     MAGEMin_C.seval("using MAGEMin_C")
-except Exception as _e:
+except Exception as _e:  # noqa: BLE001
     import warnings as _warnings
 
     _warnings.warn(
@@ -20,17 +21,43 @@
     jl = None
     jlconvert = None
 
-# Expose Public API
-from .core.base import MAGEMinBase
-from .core.engine import single_point_minimization_with_conversion
-from .core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, calculate_kd_fe_mg
+__all__ = [
+    "GarnetGenerator",
+    "MAGEMinAssemblageCalculator",
+    "MAGEMinBase",
+    "MAGEMinGarnetCalculator",
+    "MAGEMinPTGridCalculator",
+    "MAGEMin_C",
+    "PhaseFunctions",
+    "PhasePTEstimator",
+    "bulk_rock",
+    "calculate_kd_fe_mg",
+    "extract_end_member",
+    "generate_distribution",
+    "get_oxide_apfu",
+    "get_phase_chemistry",
+    "get_phase_mg_number",
+    "jl",
+    "jlconvert",
+    "phase_frac",
+    "single_point_minimization_with_conversion",
+]
 
-from .calculators.garnet import MAGEMinGarnetCalculator
+# Expose Public API
 from .calculators.assemblage import MAGEMinAssemblageCalculator
-from .calculators.pt_estimation import PhasePTEstimator
+from .calculators.garnet import MAGEMinGarnetCalculator
 from .calculators.phase_search import PhaseFunctions
+from .calculators.pt_estimation import PhasePTEstimator
 from .calculators.pt_grid import MAGEMinPTGridCalculator
-
+from .core.base import MAGEMinBase
+from .core.engine import single_point_minimization_with_conversion
+from .core.phase_properties import (
+    calculate_kd_fe_mg,
+    extract_end_member,
+    get_oxide_apfu,
+    get_phase_chemistry,
+    get_phase_mg_number,
+    phase_frac,
+)
 from .models.garnet_growth import GarnetGenerator, generate_distribution
-
 from .utils import bulk_rock
diff --git a/src/phasetools/calculators/assemblage.py b/src/phasetools/calculators/assemblage.py
index 9dc213b..d6e4134 100644
--- a/src/phasetools/calculators/assemblage.py
+++ b/src/phasetools/calculators/assemblage.py
@@ -1,8 +1,12 @@
-import numpy as np
 import sys
+
+import numpy as np
+
+from phasetools import MAGEMin_C
+
 from ..core.base import MAGEMinBase
 from ..core.phase_properties import phase_frac
-from phasetools import MAGEMin_C
+
 
 class MAGEMinAssemblageCalculator(MAGEMinBase):
     """
diff --git a/src/phasetools/calculators/garnet.py b/src/phasetools/calculators/garnet.py
index 8d2a30c..7dd76d6 100644
--- a/src/phasetools/calculators/garnet.py
+++ b/src/phasetools/calculators/garnet.py
@@ -1,10 +1,15 @@
-import numpy as np
 import sys
-from .pt_grid import MAGEMinPTGridCalculator
-from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu
-from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac
+
+import numpy as np
+from juliacall import Main as jl
+from juliacall import convert as jlconvert
+
 from phasetools import MAGEMin_C
-from juliacall import Main as jl, convert as jlconvert
+
+from ..core.phase_properties import get_oxide_apfu, phase_frac
+from ..utils.bulk_rock import atomic_frac_to_wt_frac, atomic_mass_dict
+from .pt_grid import MAGEMinPTGridCalculator
+
 
 class MAGEMinGarnetCalculator(MAGEMinPTGridCalculator):
     """High-level wrappers for garnet-focused MAGEMin calculations."""
diff --git a/src/phasetools/calculators/phase_search.py b/src/phasetools/calculators/phase_search.py
index cdfa8d9..b410ba2 100644
--- a/src/phasetools/calculators/phase_search.py
+++ b/src/phasetools/calculators/phase_search.py
@@ -1,9 +1,13 @@
-import numpy as np
 import warnings
+
+import numpy as np
 from scipy import optimize
+
+from phasetools import MAGEMin_C
+
 from ..core.base import MAGEMinBase
 from ..core.phase_properties import phase_frac
-from phasetools import MAGEMin_C
+
 
 class PhaseFunctions(MAGEMinBase):
     """Utility functions for phase stability and fractionation."""
diff --git a/src/phasetools/calculators/pt_estimation.py b/src/phasetools/calculators/pt_estimation.py
index 88227d7..edda66c 100644
--- a/src/phasetools/calculators/pt_estimation.py
+++ b/src/phasetools/calculators/pt_estimation.py
@@ -1,9 +1,12 @@
 import numpy as np
-import phasetools
 from scipy import optimize
+
+import phasetools
+
 from ..core.base import MAGEMinBase
 from ..core.phase_properties import extract_end_member, get_oxide_apfu
 
+
 class PhasePTEstimator(MAGEMinBase):
     """
     Generalized utility for estimating P-T conditions by minimizing the misfit 
@@ -57,7 +60,7 @@ def calculate_misfit(self, pt_guess, target_chemistry, phase, components, comp_t
             sigma = np.where(sigma < 1e-12, 1.0, sigma)
             
             return np.sqrt(np.mean((diff / sigma)**2))
-        except Exception:
+        except Exception:  # noqa: BLE001 — optimiser penalty: any MAGEMin failure means this P-T is infeasible
             return 1e6 
 
     def solve(self, target_chemistry, phase, components, comp_type, bounds, uncertainty=None,
@@ -117,7 +120,7 @@ def callback(x, *args):
                 sys_in=self.sys_in, rm_list=self.rm_list
             )
             res.modelled_composition = self._get_phase_composition(out, phase, components, comp_type)
-        except Exception:
+        except Exception:  # noqa: BLE001 — any calc failure → no composition to compare
             res.modelled_composition = None
 
         return res
diff --git a/src/phasetools/calculators/pt_grid.py b/src/phasetools/calculators/pt_grid.py
index 7fee7ec..cbad53f 100644
--- a/src/phasetools/calculators/pt_grid.py
+++ b/src/phasetools/calculators/pt_grid.py
@@ -1,11 +1,22 @@
-import numpy as np
 import sys
-from juliacall import Main as jl, convert as jlconvert
-from ..core.base import MAGEMinBase
-from ..core.phase_properties import extract_end_member, get_oxide_apfu, get_phase_chemistry, get_phase_mg_number, _phase_indices
-from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac
+
+import numpy as np
+from juliacall import Main as jl
+from juliacall import convert as jlconvert
+
 from phasetools import MAGEMin_C
 
+from ..core.base import MAGEMinBase
+from ..core.phase_properties import (
+    _phase_indices,
+    extract_end_member,
+    get_oxide_apfu,
+    get_phase_chemistry,
+    get_phase_mg_number,
+)
+from ..utils.bulk_rock import atomic_frac_to_wt_frac, atomic_mass_dict
+
+
 class MAGEMinPTGridCalculator(MAGEMinBase):
     """
     Calculator for performing P-T grid minimisations and extracting phase properties.
@@ -69,7 +80,7 @@ def get_all_unique_phases(self, grid_out=None):
         unique_phases = set()
         for o in out:
             unique_phases.update(o.ph)
-        return sorted(list(unique_phases))
+        return sorted(unique_phases)
 
     def get_phase_endmembers(self, phase, grid_out=None):
         """
@@ -207,6 +218,14 @@ def _precreate(prefix, names):
                 instances[k]["Fe2"] = np.full(P_len, np.nan)
                 instances[k]["Fe3"] = np.full(P_len, np.nan)
 
+        def _fill(key, value, idx):
+            vals = np.atleast_1d(np.asarray(value, dtype=float))
+            n = len(vals)
+            for k in range(n_inst):
+                if k < n:
+                    instances[k][key][idx] = vals[k]
+                # else stays NaN
+
         for i in range(P_len):
             if phase not in out[i].ph:
                 continue
@@ -217,32 +236,24 @@ def _precreate(prefix, names):
                 instances[k]["vol_frac"][i] = float(out[i].ph_frac_vol[j])
             # instances k >= len(n_idx) stay NaN
 
-            def _fill(key, value):
-                vals = np.atleast_1d(np.asarray(value, dtype=float))
-                n = len(vals)
-                for k in range(n_inst):
-                    if k < n:
-                        instances[k][key][i] = vals[k]
-                    # else stays NaN
-
             if end_members:
                 for em in end_members:
-                    _fill(f"em_{em}", extract_end_member(phase, out[i], em, self.sys_in, instance='all'))
+                    _fill(f"em_{em}", extract_end_member(phase, out[i], em, self.sys_in, instance='all'), i)
             if oxides:
                 apfu = get_oxide_apfu(out[i], phase, oxides, instance='all')
-                for ox in oxides: _fill(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0)))
+                for ox in oxides: _fill(f"ox_apfu_{ox}", apfu.get(ox, np.zeros(0)), i)
             if chemistry:
                 chem = get_phase_chemistry(out[i], phase, chemistry, self.sys_in, instance='all')
-                for ox in chemistry: _fill(f"chem_{ox}", chem.get(ox, np.zeros(0)))
+                for ox in chemistry: _fill(f"chem_{ox}", chem.get(ox, np.zeros(0)), i)
             if cations:
                 cat_vals = self._extract_cations_from_apfu(out[i], phase, cations, self.sys_in)
-                for c in cations: _fill(f"cat_{c}", cat_vals[f"cat_{c}"])
+                for c in cations: _fill(f"cat_{c}", cat_vals[f"cat_{c}"], i)
             if mg_number:
-                _fill("Mg_number", get_phase_mg_number(out[i], phase, instance='all'))
+                _fill("Mg_number", get_phase_mg_number(out[i], phase, instance='all'), i)
             if fe_split:
                 split = self._extract_fe_split_from_apfu(out[i], phase, instance='all')
-                _fill("Fe2", split["Fe2"])
-                _fill("Fe3", split["Fe3"])
+                _fill("Fe2", split["Fe2"], i)
+                _fill("Fe3", split["Fe3"], i)
 
         return instances
 
diff --git a/src/phasetools/core/base.py b/src/phasetools/core/base.py
index d81bd75..cec49b2 100644
--- a/src/phasetools/core/base.py
+++ b/src/phasetools/core/base.py
@@ -1,10 +1,12 @@
-import numpy as np
-from juliacall import Main as jl, convert as jlconvert
+from juliacall import Main as jl
+from juliacall import convert as jlconvert
 from molmass import Formula
+
 from phasetools import MAGEMin_C
 
 from ..utils.bulk_rock import convert_mol_percent_to_wt_percent, get_molar_mass_dict
 
+
 class MAGEMinBase:
     def __init__(self, db="ig", dataset=636, verbose=False):
         self.db = db
@@ -27,7 +29,6 @@ def _copy_state_to(self, other):
         other._stoich_map = self._stoich_map
 
     def setup_bulk_composition(self, Xoxides, X, sys_in, rm_list=None):
-        from ..utils.bulk_rock import convert_mol_percent_to_wt_percent, get_molar_mass_dict
 
         # Store original inputs for later modification (e.g. in estimators)
         self.original_Xoxides = list(Xoxides)
@@ -72,7 +73,7 @@ def setup_bulk_composition(self, Xoxides, X, sys_in, rm_list=None):
                         if ox not in self._stoich_map:
                             self._stoich_map[ox] = {}
                         self._stoich_map[ox][el] = float(df.loc[el, 'Count'])
-            except:
+            except (ValueError, KeyError, TypeError, AttributeError):
                 pass
 
     def _extract_fe_split_from_apfu(self, out, phase, instance=0):
diff --git a/src/phasetools/core/engine.py b/src/phasetools/core/engine.py
index d471730..6f979b9 100644
--- a/src/phasetools/core/engine.py
+++ b/src/phasetools/core/engine.py
@@ -1,6 +1,9 @@
-from juliacall import Main as jl, convert as jlconvert
+from juliacall import Main as jl
+from juliacall import convert as jlconvert
+
 from phasetools import MAGEMin_C
 
+
 def single_point_minimization_with_conversion(P, T, data, X, Xoxides, sys_in="wt", rm_list=None):
     """
     Perform single-point minimisation with automatic Python-to-Julia conversion, including optional removal list.
diff --git a/src/phasetools/core/phase_properties.py b/src/phasetools/core/phase_properties.py
index 4a78141..b938653 100644
--- a/src/phasetools/core/phase_properties.py
+++ b/src/phasetools/core/phase_properties.py
@@ -1,6 +1,7 @@
-import numpy as np
 import warnings
 
+import numpy as np
+
 
 def _phase_indices(out: object, phase: str, instance: int | str = 0) -> list[int]:
     """Return the index/indices of ``phase`` in ``out.ph``.
@@ -29,7 +30,7 @@ def _phase_indices(out: object, phase: str, instance: int | str = 0) -> list[int
     if instance == 'all':
         return idx
     if not isinstance(instance, (int, np.integer)):
-        raise ValueError(f"instance must be an integer index or 'all', got {instance!r}")
+        raise TypeError(f"instance must be an integer index or 'all', got {instance!r}")
     if not idx:
         return []
     if instance < 0:
@@ -178,7 +179,7 @@ def phase_frac(phase, MAGEMinOutput, sys_in):
         if not found:
             return 0.0
         return float(total)
-    except:
+    except (ValueError, IndexError, AttributeError, TypeError):
         return 0.0
 
 def get_phase_mg_number(out, ph, instance=0):
@@ -281,7 +282,7 @@ def get_phase_mg2_number(out: object, ph: str, instance: int | str = 0) -> float
         if instance != 'all':
             return float(mg / denominator) if denominator > 0 else 0.0
         return np.where(denominator > 0, mg / denominator, 0.0)
-    except Exception:
+    except (ValueError, IndexError, AttributeError, KeyError, TypeError):
         return 0.0
 
 def calculate_kd_fe_mg(out, phase1, phase2, use_fe2_only=False):
diff --git a/src/phasetools/julia_setup.py b/src/phasetools/julia_setup.py
index e4999f9..c40d144 100644
--- a/src/phasetools/julia_setup.py
+++ b/src/phasetools/julia_setup.py
@@ -5,7 +5,7 @@
 
 
 def _run_julia(args):
-    return subprocess.run(["julia", *args], capture_output=True, text=True)
+    return subprocess.run(["julia", *args], capture_output=True, text=True, check=False)
 
 
 def check_julia_installed():
@@ -68,7 +68,7 @@ def main(argv=None):
         try:
             install_magemin()
             print("MAGEMin_C installed successfully.")
-        except Exception as exc:
+        except (RuntimeError, OSError) as exc:
             print(f"Error: {exc}", file=sys.stderr)
             return 1
 
diff --git a/src/phasetools/models/garnet_growth.py b/src/phasetools/models/garnet_growth.py
index bad72e9..f5e160b 100644
--- a/src/phasetools/models/garnet_growth.py
+++ b/src/phasetools/models/garnet_growth.py
@@ -1,11 +1,15 @@
-import numpy as np
+from typing import Any
+
 import matplotlib.pyplot as plt
+import numpy as np
+from juliacall import Main as jl
+from juliacall import convert as jlconvert
 from scipy.interpolate import interp1d
 from scipy.stats import norm
-from typing import Any
-from juliacall import Main as jl, convert as jlconvert
+
 from ..calculators.garnet import MAGEMinGarnetCalculator
 
+
 def generate_distribution(n_classes, r_min, dr, fnr, Gn, tGn):
     """
     Generates a distribution of radial sizes and associated garnet volumes and formation times.
@@ -181,7 +185,7 @@ def _get_size_distribution(self, size_dist, r):
                 raise ValueError("User-defined distribution must have length equal to garnet_classes")
             finp = user_dist
         else:
-            raise ValueError("size_dist must be a string ('N' or 'U') or a numeric array")
+            raise TypeError("size_dist must be a string ('N' or 'U') or a numeric array")
         return finp
 
     def _normalise_distribution(self, finp, r):
@@ -415,15 +419,15 @@ def generate_garnets(self, size_dist='N'):
 
         GVG = GVn[ind]
 
-        tG, TG, PG, MnG, MgG, FeG, CaG = self._slice_arrays(ind)
+        tG, TG, PG, MnG, MgG, FeG, _CaG = self._slice_arrays(ind)
 
         # Generate radius classes and distributions
-        n_classes, r, dr, finp, fnr = self._build_size_distribution(size_dist)
+        n_classes, _r, dr, _finp, fnr = self._build_size_distribution(size_dist)
 
         Gn = GVG / np.max(GVG)
         tGn = tG
 
-        G, t_arr, r_r, R = generate_distribution(n_classes, self.r_min, dr, fnr, Gn, tGn)
+        _G, t_arr, _r_r, R = generate_distribution(n_classes, self.r_min, dr, fnr, Gn, tGn)
 
         # Interpolate physical properties along the garnet growth
         PGrw = self._interp(tG, PG, t_arr)
@@ -431,7 +435,6 @@ def generate_garnets(self, size_dist='N'):
         Mnrw = self._interp(tG, MnG, t_arr)
         Mgrw = self._interp(tG, MgG, t_arr)
         Ferw = self._interp(tG, FeG, t_arr)
-        Carw = 1 - Mnrw - Mgrw - Ferw
 
         garnets = []
         for i in range(n_classes):
@@ -443,7 +446,6 @@ def generate_garnets(self, size_dist='N'):
             Mnr1 = Mnrw[ind_range]
             Mgr1 = Mgrw[ind_range]
             Fer1 = Ferw[ind_range]
-            Car1 = 1 - Mnr1 - Mgr1 - Fer1
 
             dRr = Rr1[-1] / self.nR_diff
             Rrz = np.arange(dRr, Rr1[-1] + dRr, dRr)
@@ -501,7 +503,7 @@ def plot_garnet_summary(self, size_dist='N', garnet_no=0, path=None, plot_fig=Tr
 
         GVG = GVn[ind]
 
-        tG, TG, PG, MnG, MgG, FeG, CaG = self._slice_arrays(ind)
+        tG, TG, PG, MnG, MgG, FeG, _CaG = self._slice_arrays(ind)
 
         n_classes, r, dr, finp, fnr = self._build_size_distribution(size_dist)
 
@@ -509,7 +511,7 @@ def plot_garnet_summary(self, size_dist='N', garnet_no=0, path=None, plot_fig=Tr
         tGn = tG
 
         G, t_arr, r_r, R = generate_distribution(n_classes, self.r_min, dr, fnr, Gn, tGn)
-        
+
         # Interpolate physical properties along garnet growth
         PGrw = self._interp(tG, PG, t_arr)
         TGrw = self._interp(tG, TG, t_arr)
@@ -621,7 +623,7 @@ def plot_garnet_summary(self, size_dist='N', garnet_no=0, path=None, plot_fig=Tr
         except IndexError:
             last_zero_idx = -1
         ind = np.arange(last_zero_idx+1, first_one_idx+1)
-        tG, TG, PG, MnG, MgG, FeG, CaG = self._slice_arrays(ind)
+        tG, TG, PG, MnG, MgG, FeG, _CaG = self._slice_arrays(ind)
         n_classes, r, dr, finp, fnr = self._build_size_distribution(size_dist)
         Gn = GVn[ind] / np.max(GVn[ind])
         G, t_arr, r_r, R = generate_distribution(n_classes, self.r_min, dr, fnr, Gn, tG)
diff --git a/src/phasetools/models/magma_ocean.py b/src/phasetools/models/magma_ocean.py
index 4106300..2ccb430 100644
--- a/src/phasetools/models/magma_ocean.py
+++ b/src/phasetools/models/magma_ocean.py
@@ -1,12 +1,16 @@
-import numpy as np
-import sys
 import warnings
+from typing import Any
+
+import numpy as np
+from juliacall import Main as jl
+from juliacall import convert as jlconvert
 from scipy import optimize
-from typing import List, Dict, Any, Tuple, Optional
-from ..core.base import MAGEMinBase
-from ..core.phase_properties import phase_frac, get_phase_chemistry
+
 from phasetools import MAGEMin_C
-from juliacall import Main as jl, convert as jlconvert
+
+from ..core.base import MAGEMinBase
+from ..core.phase_properties import phase_frac
+
 
 class MagmaOcean(MAGEMinBase):
     """
@@ -75,8 +79,10 @@ def get_volume_between_radii(self, r1: float, r2: float) -> float:
         """Calculate volume of a spherical shell between radii r1 and r2."""
         return (4.0/3.0) * np.pi * (np.abs(r1**3 - r2**3))
 
-    def find_temperature_at_vol_frac(self, P: float, target_vol_frac: float, bracket: List[float] = [800, 3000]) -> float:
+    def find_temperature_at_vol_frac(self, P: float, target_vol_frac: float, bracket: list[float] | None = None) -> float:
         """Find the temperature at which a specific volume fraction of solid is reached."""
+        if bracket is None:
+            bracket = [800, 3000]
         def func(T):
             out = MAGEMin_C.single_point_minimization(P, T, self.data, X=self.X, Xoxides=self.Xoxides, sys_in=self.sys_in, rm_list=self.rm_list)
             # Total solid fraction = 1 - melt fraction (liq)
@@ -112,7 +118,7 @@ def get_phase_chemistry_at_index(self, out, i: int) -> np.ndarray:
         else:
             return np.array(obj.Comp, dtype=float)
 
-    def run_stage_0(self, p_start: float, p_end: float, solid_frac: float = 0.5, p_intervals: int = 20) -> Tuple[Dict[str, Any], np.ndarray]:
+    def run_stage_0(self, p_start: float, p_end: float, solid_frac: float = 0.5, p_intervals: int = 20) -> tuple[dict[str, Any], np.ndarray]:
         """
         Stage 0: Equilibrium crystallisation from 0 to target solid fraction.
         """
@@ -171,14 +177,14 @@ def run_stage_0(self, p_start: float, p_end: float, solid_frac: float = 0.5, p_i
         
         return results, avg_melt
 
-    def run_fractional_stages(self, 
-                               starting_melt: np.ndarray, 
-                               p_start: float, 
-                               p_end: float, 
+    def run_fractional_stages(self,
+                               starting_melt: np.ndarray,
+                               p_start: float,
+                               p_end: float,
                                vol_step: float = 0.05,
                                starting_vol_frac: float = 0.5,
                                n_stages: int = 10,
-                               float_phases: List[str] = ['pl', 'q', 'fsp', 'san', 'ksp']) -> List[Dict[str, Any]]:
+                               float_phases: list[str] | None = None) -> list[dict[str, Any]]:
         """
         Fractional crystallisation stages of the remaining melt.
         
@@ -203,6 +209,8 @@ def run_fractional_stages(self,
         float_phases : list[str]
             List of phases that float to form the crust.
         """
+        if float_phases is None:
+            float_phases = ['pl', 'q', 'fsp', 'san', 'ksp']
         all_stage_results = []
         current_melt_comp = starting_melt
 
diff --git a/tests/test_bulk_rock.py b/tests/test_bulk_rock.py
index 83cc5da..e8db5e0 100644
--- a/tests/test_bulk_rock.py
+++ b/tests/test_bulk_rock.py
@@ -1,13 +1,15 @@
 import unittest
+
 import numpy as np
 
 from phasetools.utils.bulk_rock import (
-    get_molar_mass_dict,
     get_atomic_mass_dict,
+    get_molar_mass_dict,
     mol_fractions_to_wt_fractions,
     wt_fractions_to_mol_fractions,
 )
 
+
 class TestFractionConverters(unittest.TestCase):
     """Tests for the non-normalising mole/weight fraction converters."""
 
diff --git a/tests/test_garnet.py b/tests/test_garnet.py
index 326af8d..5bb0828 100644
--- a/tests/test_garnet.py
+++ b/tests/test_garnet.py
@@ -11,9 +11,10 @@
 """
 
 import unittest
-import numpy as np
 from unittest.mock import MagicMock, patch
 
+import numpy as np
+
 from phasetools.calculators.garnet import MAGEMinGarnetCalculator
 
 
@@ -66,10 +67,12 @@ def test_absent_garnet_returns_zeros(self):
 
     def test_invalid_basis_raises(self):
         """Unsupported fe_basis values must raise ValueError at construction."""
-        with patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__',
-                   return_value=None):
-            with self.assertRaises(ValueError):
-                MAGEMinGarnetCalculator(db='ig', fe_basis='Fe3')
+        with (
+            patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__',
+                  return_value=None),
+            self.assertRaises(ValueError),
+        ):
+            MAGEMinGarnetCalculator(db='ig', fe_basis='Fe3')
 
     def test_default_is_feot(self):
         """The default fe_basis is 'FeOt' (community convention)."""
diff --git a/tests/test_grid_extraction.py b/tests/test_grid_extraction.py
index ac35947..e6b8e7f 100644
--- a/tests/test_grid_extraction.py
+++ b/tests/test_grid_extraction.py
@@ -20,14 +20,19 @@
 
 import unittest
 import warnings
-import numpy as np
 from unittest.mock import MagicMock, patch
 
+import numpy as np
+
 from phasetools.calculators.pt_grid import MAGEMinPTGridCalculator
 from phasetools.core.phase_properties import (
-    get_oxide_apfu, get_phase_chemistry, extract_end_member,
-    get_phase_mg_number, get_phase_mg2_number, get_phase_fe_split,
     _phase_indices,
+    extract_end_member,
+    get_oxide_apfu,
+    get_phase_chemistry,
+    get_phase_fe_split,
+    get_phase_mg2_number,
+    get_phase_mg_number,
 )
 
 # mpe-style oxide ordering
@@ -105,7 +110,7 @@ def test_absent_phase_empty(self):
         self.assertEqual(_phase_indices(self.out, 'ep', 0), [])
 
     def test_invalid_instance_type_raises(self):
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             _phase_indices(self.out, 'dio', 'first')
 
 
diff --git a/tests/test_magma_ocean.py b/tests/test_magma_ocean.py
index 3dc0143..0adb3f5 100644
--- a/tests/test_magma_ocean.py
+++ b/tests/test_magma_ocean.py
@@ -1,6 +1,8 @@
 import unittest
+
 from phasetools.models.magma_ocean import MagmaOcean
 
+
 class TestLMOFix(unittest.TestCase):
     """
     Verification of the Lunar Magma Ocean (LMO) crystallisation logic.
@@ -20,7 +22,7 @@ def test_full_crystallisation_path(self):
         """Test the 10-stage fractional crystallisation path."""
         # Stage 0: Equilibrium (0-50 vol.%)
         # Using p_start=45.0 (core) to p_end=17.3 (approx. 50% volume radius)
-        results_0, melt_0 = self.mo.run_stage_0(p_start=45.0, p_end=17.3, solid_frac=0.5, p_intervals=5)
+        _, melt_0 = self.mo.run_stage_0(p_start=45.0, p_end=17.3, solid_frac=0.5, p_intervals=5)
 
         # Fractional Stages (Stages 1-10, 5% total LMO each)
         results_frac = self.mo.run_fractional_stages(melt_0, p_start=17.3, p_end=0.01, vol_step=0.05, n_stages=10)
diff --git a/tests/test_phase_properties.py b/tests/test_phase_properties.py
index b671a4a..be24996 100644
--- a/tests/test_phase_properties.py
+++ b/tests/test_phase_properties.py
@@ -1,14 +1,14 @@
+# Mock MAGEMin_C before importing PhaseTools components that depend on it
+import sys
 import unittest
-import numpy as np
 from unittest.mock import MagicMock, patch
 
-# Mock MAGEMin_C before importing PhaseTools components that depend on it
-import sys
 mock_magemin = MagicMock()
 sys.modules['phasetools.MAGEMin_C'] = mock_magemin
 
 from phasetools.core.base import MAGEMinBase
 
+
 class TestRedoxLogic(unittest.TestCase):
     def setUp(self):
         # Prevent Initialize_MAGEMin from calling Julia
diff --git a/tests/test_phase_properties_live.py b/tests/test_phase_properties_live.py
index 1dac3af..f133e12 100644
--- a/tests/test_phase_properties_live.py
+++ b/tests/test_phase_properties_live.py
@@ -1,10 +1,17 @@
 import unittest
+from typing import ClassVar
+
 import numpy as np
+
 from phasetools.calculators.pt_grid import MAGEMinPTGridCalculator
-from phasetools.core.phase_properties import get_phase_fe_split, get_phase_mg_number, get_phase_mg2_number, phase_frac
+from phasetools.core.phase_properties import (
+    get_phase_fe_split,
+    get_phase_mg2_number,
+    phase_frac,
+)
 
 try:
-    from juliacall import Main as jl
+    import juliacall  # noqa: F401 — side-effect import to test availability
     HAS_JULIA = True
 except ImportError:
     HAS_JULIA = False
@@ -24,7 +31,7 @@ def setUpClass(cls):
         # 30 kbar (3 GPa), 700 C - where g and dio are stable in mpe
         cls.out = cls.calc.calculate_grid(30.0, 700.0)[0]
 
-    PT_POINTS = [
+    PT_POINTS: ClassVar[list[tuple[float, float]]] = [
         (5.0, 450.0),
         (10.0, 550.0),
         (20.0, 650.0),
@@ -119,8 +126,7 @@ def test_amphibole_heuristic(self):
         fe2_h, fe3_h = split['Fe2'], split['Fe3']
         
         ph_idx = out.ph.index('amp')
-        em = {str(n): float(f) for n, f in zip(out.SS_vec[ph_idx].emNames, out.SS_vec[ph_idx].emFrac)}
-        
+
         # Sum Fe2+ and Fe3+ from endmembers
         # Use sum constraint as primary check: Fe2+ + Fe3+ = total Fe
         ox_names = [str(o) for o in out.oxides]
@@ -265,7 +271,7 @@ def setUpClass(cls):
         cls.calc = MAGEMinPTGridCalculator(db=cls.db, dataset=636)
         cls.calc.setup_bulk_composition(cls.Xoxides, cls.X, sys_in='wt')
 
-    PT_POINTS = [
+    PT_POINTS: ClassVar[list[tuple[float, float]]] = [
         (15.0, 600.0),
         (20.0, 700.0),
         (25.0, 750.0),
@@ -300,7 +306,10 @@ def test_heuristic_mg_matches_site_fractions(self):
 
     def test_kd_mg_independent_of_units(self):
         """K_D and Mg# are identical whether bulk is wt% or mol%."""
-        from phasetools.utils.bulk_rock import convert_wt_percent_to_mol_percent, get_molar_mass_dict
+        from phasetools.utils.bulk_rock import (
+            convert_wt_percent_to_mol_percent,
+            get_molar_mass_dict,
+        )
         mass_dict = get_molar_mass_dict()
         X_mol = convert_wt_percent_to_mol_percent(self.X, self.Xoxides, mass_dict)
 
@@ -355,7 +364,7 @@ def setUpClass(cls):
         cls.calc = MAGEMinPTGridCalculator(db=cls.db, dataset=636)
         cls.calc.setup_bulk_composition(cls.Xoxides, cls.X, sys_in='wt')
 
-    PT_POINTS = [
+    PT_POINTS: ClassVar[list[tuple[float, float]]] = [
         (20.0, 700.0),
         (25.0, 750.0),
         (30.0, 800.0),
diff --git a/tests/test_phase_search.py b/tests/test_phase_search.py
index 986a554..4f52447 100644
--- a/tests/test_phase_search.py
+++ b/tests/test_phase_search.py
@@ -9,9 +9,9 @@
 """
 
 import unittest
-import numpy as np
-from unittest.mock import MagicMock, patch, PropertyMock
+from unittest.mock import MagicMock
 
+import numpy as np
 
 # ---------------------------------------------------------------------------
 # Helpers