From 1ec01a172e1a3970486cc5a91c98d58b354973ee Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@cphyc.me>
Date: Fri, 17 May 2024 14:29:58 +0100
Subject: [PATCH 01/13] Simplify handling of conformal birth times in RAMSES

Also read from 'birthXXXX.outYYYYY' files if found
---
 yt/frontends/ramses/data_structures.py   |  9 ++-
 yt/frontends/ramses/fields.py            |  2 +-
 yt/frontends/ramses/io.py                | 62 ++++++---------
 yt/frontends/ramses/particle_handlers.py | 96 ++++++++++++++++++++++++
 4 files changed, 123 insertions(+), 46 deletions(-)

diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py
index 5b811cf07b2..51206aa68b8 100644
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -640,10 +640,6 @@ def _detect_output_fields(self):
             dsl.update(set(ph.field_offsets.keys()))
 
         self.particle_field_list = list(dsl)
-        cosmo = self.ds.cosmological_simulation
-        for f in self.particle_field_list[:]:
-            if f[1] == "particle_birth_time" and cosmo:
-                self.particle_field_list.append((f[0], "conformal_birth_time"))
 
         # Get the detected fields
         dsl = set()
@@ -809,6 +805,7 @@ def __init__(
         max_level=None,
         max_level_convention=None,
         default_species_fields=None,
+        use_conformal_time=None,
     ):
         # Here we want to initiate a traceback, if the reader is not built.
         if isinstance(fields, str):
@@ -833,6 +830,10 @@ def __init__(
         self._extra_particle_fields = extra_particle_fields
         self.force_cosmological = cosmological
         self._bbox = bbox
+        if use_conformal_time is not None:
+            self.use_conformal_time = use_conformal_time
+        else:
+            self.use_conformal_time = cosmological
 
         self._force_max_level = self._sanitize_max_level(
             max_level, max_level_convention
diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py
index 4c871ec08b8..989d112860a 100644
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -180,7 +180,7 @@ def setup_particle_fields(self, ptype):
         super().setup_particle_fields(ptype)
 
         def star_age(field, data):
-            if data.ds.cosmological_simulation:
+            if data.ds.cosmological_simulation and data.ds.use_conformal_time:
                 conformal_age = data[ptype, "conformal_birth_time"]
                 physical_age = convert_ramses_conformal_time_to_physical_age(
                     data.ds, conformal_age
diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index 546ba333f5c..0afa1464a9e 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -75,7 +75,7 @@ def convert_ramses_conformal_time_to_physical_age(
     return (tsim - t) * t_scale
 
 
-def _ramses_particle_binary_file_handler(particle, subset, fields, count):
+def _ramses_particle_binary_file_handler(particle_handler, subset, fields, count):
     """General file handler for binary file, called by _read_particle_subset
 
     Parameters
@@ -91,10 +91,9 @@ def _ramses_particle_binary_file_handler(particle, subset, fields, count):
     """
     tr = {}
     ds = subset.domain.ds
-    current_time = ds.current_time.in_units("code_time").v
-    foffsets = particle.field_offsets
-    fname = particle.fname
-    data_types = particle.field_types
+    foffsets = particle_handler.field_offsets
+    fname = particle_handler.fname
+    data_types = particle_handler.field_types
     with FortranFile(fname) as fd:
         # We do *all* conversion into boxlen here.
         # This means that no other conversions need to be applied to convert
@@ -103,24 +102,23 @@ def _ramses_particle_binary_file_handler(particle, subset, fields, count):
             if count == 0:
                 tr[field] = np.empty(0, dtype=data_types[field])
                 continue
-            fd.seek(foffsets[field])
-            dt = data_types[field]
-            tr[field] = fd.read_vector(dt)
+            # Sentinel value: -1 means we don't have this field
+            if foffsets[field] == -1:
+                tr[field] = np.empty(count, dtype=data_types[field])
+            else:
+                fd.seek(foffsets[field])
+                dt = data_types[field]
+                tr[field] = fd.read_vector(dt)
             if field[1].startswith("particle_position"):
                 np.divide(tr[field], ds["boxlen"], tr[field])
-            if ds.cosmological_simulation and field[1] == "particle_birth_time":
-                conformal_time = tr[field]
-                physical_age = convert_ramses_conformal_time_to_physical_age(
-                    ds, conformal_time
-                )
-                tr[field] = current_time - physical_age
-                # arbitrarily set particles with zero conformal_age to zero
-                # particle_age. This corresponds to DM particles.
-                tr[field][conformal_time == 0] = 0
+
+            # Handle over to field handler for special cases, like particle_birth_times
+            particle_handler.handle_field(field, tr)
+
     return tr
 
 
-def _ramses_particle_csv_file_handler(particle, subset, fields, count):
+def _ramses_particle_csv_file_handler(particle_handler, subset, fields, count):
     """General file handler for csv file, called by _read_particle_subset
 
     Parameters
@@ -138,9 +136,8 @@ def _ramses_particle_csv_file_handler(particle, subset, fields, count):
 
     tr = {}
     ds = subset.domain.ds
-    current_time = ds.current_time.in_units("code_time").v
-    foffsets = particle.field_offsets
-    fname = particle.fname
+    foffsets = particle_handler.field_offsets
+    fname = particle_handler.fname
 
     list_field_ind = [
         (field, foffsets[field]) for field in sorted(fields, key=lambda a: foffsets[a])
@@ -159,15 +156,9 @@ def _ramses_particle_csv_file_handler(particle, subset, fields, count):
         tr[field] = dat[ind].to_numpy()
         if field[1].startswith("particle_position"):
             np.divide(tr[field], ds["boxlen"], tr[field])
-        if ds.cosmological_simulation and field[1] == "particle_birth_time":
-            conformal_time = tr[field]
-            physical_age = convert_ramses_conformal_time_to_physical_age(
-                ds, conformal_time
-            )
-            tr[field] = current_time - physical_age
-            # arbitrarily set particles with zero conformal_age to zero
-            # particle_age. This corresponds to DM particles.
-            tr[field][conformal_time == 0] = 0
+
+        particle_handler.handle_field(field, tr)
+
     return tr
 
 
@@ -294,23 +285,12 @@ def _read_particle_subset(self, subset, fields):
             ok = False
             for ph in subset.domain.particle_handlers:
                 if ph.ptype == ptype:
-                    foffsets = ph.field_offsets
-                    data_types = ph.field_types
                     ok = True
                     count = ph.local_particle_count
                     break
             if not ok:
                 raise YTFieldTypeNotFound(ptype)
 
-            cosmo = self.ds.cosmological_simulation
-            if (ptype, "particle_birth_time") in foffsets and cosmo:
-                foffsets[ptype, "conformal_birth_time"] = foffsets[
-                    ptype, "particle_birth_time"
-                ]
-                data_types[ptype, "conformal_birth_time"] = data_types[
-                    ptype, "particle_birth_time"
-                ]
-
             tr.update(ph.reader(subset, subs_fields, count))
 
         return tr
diff --git a/yt/frontends/ramses/particle_handlers.py b/yt/frontends/ramses/particle_handlers.py
index 2b42186583c..32da4d1aca3 100644
--- a/yt/frontends/ramses/particle_handlers.py
+++ b/yt/frontends/ramses/particle_handlers.py
@@ -2,6 +2,7 @@
 import os
 import struct
 from collections.abc import Callable
+from functools import cached_property
 from itertools import chain, count
 from typing import TYPE_CHECKING, Any
 
@@ -18,6 +19,7 @@
     _ramses_particle_csv_file_handler,
     _read_part_binary_file_descriptor,
     _read_part_csv_file_descriptor,
+    convert_ramses_conformal_time_to_physical_age,
 )
 
 if TYPE_CHECKING:
@@ -156,6 +158,26 @@ def header(self) -> dict[str, Any]:
         self.read_header()
         return self._header
 
+    def handle_field(
+        self, field: tuple[str, str], data_dict: dict[tuple[str, str], np.ndarray]
+    ):
+        """
+        This function allows custom code to be called to handle special cases,
+        such as the particle birth time.
+
+        It updates the `data_dict` dictionary with the new data.
+
+        Parameters
+        ----------
+        field : tuple[str, str]
+            The field name.
+        data_dict : dict[tuple[str, str], np.ndarray]
+            A dictionary containing the data.
+
+        By default, this function does nothing.
+        """
+        pass
+
 
 _default_dtypes: dict[int, str] = {
     1: "c",  # char
@@ -293,9 +315,67 @@ def build_iterator():
             )
             self.ds._warned_extra_fields["io"] = True
 
+        if self.has_birth_file:
+            # Sentinel value: this will prevent birth times from being read from
+            # the partXXXXX.outYYYYY files
+            field_offsets[ptype, "particle_birth_time"] = -1
+            _pfields[ptype, "particle_birth_time"] = "d"
+            self.ds.use_conformal_time = False
+        elif (
+            self.ds.use_conformal_time
+            and (ptype, "particle_birth_time") in field_offsets
+        ):
+            field_offsets[ptype, "conformal_birth_time"] = field_offsets[
+                ptype, "particle_birth_time"
+            ]
+            _pfields[ptype, "conformal_birth_time"] = _pfields[
+                ptype, "particle_birth_time"
+            ]
+
         self._field_offsets = field_offsets
         self._field_types = _pfields
 
+    @property
+    def birth_file_fname(self):
+        basename = os.path.abspath(self.ds.root_folder)
+        iout = int(
+            os.path.basename(self.ds.parameter_filename).split(".")[0].split("_")[1]
+        )
+        icpu = self.domain_id
+
+        fname = os.path.join(basename, f"birth_{iout:05d}.out{icpu:05d}")
+        return fname
+
+    @cached_property
+    def has_birth_file(self):
+        return os.path.exists(self.birth_file_fname)
+
+    def handle_field(
+        self, field: tuple[str, str], data_dict: dict[tuple[str, str], np.ndarray]
+    ):
+        _ptype, fname = field
+        if not (fname == "particle_birth_time" and self.ds.cosmological_simulation):
+            return
+
+        # If the birth files exist, read from them
+        if self.has_birth_file:
+            with FortranFile(self.birth_file_fname) as fd:
+                # Note: values are written in Gyr, so we need to convert back to code_time
+                data_dict[field] = (
+                    self.ds.arr(fd.read_vector("d"), "Gyr").to("code_time").v
+                )
+
+            return
+
+        # Otherwise, convert conformal time to physical age
+        ds = self.ds
+        conformal_time = data_dict[field]
+        physical_age = convert_ramses_conformal_time_to_physical_age(ds, conformal_time)
+        current_time = ds.current_time.in_units("code_time").v
+        # arbitrarily set particles with zero conformal_age to zero
+        # particle_age. This corresponds to DM particles.
+        data_dict[field] = np.where(conformal_time > 0, current_time - physical_age, 0)
+
 
 class SinkParticleFileHandler(ParticleFileHandler):
     """Handle sink files"""
@@ -414,3 +494,19 @@ def read_header(self):
 
         self._field_offsets = field_offsets
         self._field_types = _pfields
+
+    def handle_field(
+        self, field: tuple[str, str], data_dict: dict[tuple[str, str], np.ndarray]
+    ):
+        _ptype, fname = field
+        if not (fname == "particle_birth_time" and self.ds.cosmological_simulation):
+            return
+
+        # onvert conformal time to physical age
+        ds = self.ds
+        conformal_time = data_dict[field]
+        physical_age = convert_ramses_conformal_time_to_physical_age(ds, conformal_time)
+        current_time = ds.current_time.in_units("code_time").v
+        # arbitrarily set particles with zero conformal_age to zero
+        # particle_age. This corresponds to DM particles.
+        data_dict[field] = np.where(conformal_time > 0, current_time - physical_age, 0)

From a4e0068c66544208c82ba362282e1af0a58578ff Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@cphyc.me>
Date: Thu, 2 Nov 2023 17:04:40 +0000
Subject: [PATCH 02/13] Use analytical integration rather than approximate
 Euler method

---
 yt/frontends/ramses/data_structures.py   |  58 ++++-------
 yt/frontends/ramses/fields.py            |  39 +++++--
 yt/frontends/ramses/io.py                |   6 ++
 yt/frontends/ramses/particle_handlers.py |  14 +--
 yt/utilities/lib/cosmology_time.pyx      | 127 ++++++++++++++++++++---
 5 files changed, 172 insertions(+), 72 deletions(-)

diff --git a/yt/frontends/ramses/data_structures.py b/yt/frontends/ramses/data_structures.py
index 51206aa68b8..6658ee91c07 100644
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -16,10 +16,9 @@
 from yt.geometry.oct_container import RAMSESOctreeContainer
 from yt.geometry.oct_geometry_handler import OctreeIndex
 from yt.utilities.cython_fortran_utils import FortranFile as fpu
-from yt.utilities.lib.cosmology_time import friedman
+from yt.utilities.lib.cosmology_time import t_frw, tau_frw
 from yt.utilities.on_demand_imports import _f90nml as f90nml
 from yt.utilities.physical_constants import kb, mp
-from yt.utilities.physical_ratios import cm_per_mpc
 
 from .definitions import (
     OUTPUT_DIR_EXP,
@@ -830,10 +829,6 @@ def __init__(
         self._extra_particle_fields = extra_particle_fields
         self.force_cosmological = cosmological
         self._bbox = bbox
-        if use_conformal_time is not None:
-            self.use_conformal_time = use_conformal_time
-        else:
-            self.use_conformal_time = cosmological
 
         self._force_max_level = self._sanitize_max_level(
             max_level, max_level_convention
@@ -880,6 +875,16 @@ def __init__(
             if FH.any_exist(self):
                 self.fluid_types += (FH.ftype,)
 
+        if use_conformal_time is not None:
+            self.use_conformal_time = use_conformal_time
+        elif self.cosmological_simulation:
+            if "rt" in self.fluid_types:
+                self.use_conformal_time = False
+            else:
+                self.use_conformal_time = True
+        else:
+            self.use_conformal_time = False
+
         self.storage_filename = storage_filename
 
     @staticmethod
@@ -1056,15 +1061,14 @@ def caster(val):
             is_cosmological = not (
                 rheader["time"] >= 0 and rheader["H0"] == 1 and rheader["aexp"] == 1
             )
-
         if not is_cosmological:
-            self.cosmological_simulation = 0
+            self.cosmological_simulation = False
             self.current_redshift = 0
             self.hubble_constant = 0
             self.omega_matter = 0
             self.omega_lambda = 0
         else:
-            self.cosmological_simulation = 1
+            self.cosmological_simulation = True
             self.current_redshift = (1.0 / rheader["aexp"]) - 1.0
             self.omega_lambda = rheader["omega_l"]
             self.omega_matter = rheader["omega_m"]
@@ -1075,38 +1079,14 @@ def caster(val):
             force_max_level += self.min_level + 1
         self.max_level = min(force_max_level, rheader["levelmax"]) - self.min_level - 1
 
-        if self.cosmological_simulation == 0:
+        if not self.cosmological_simulation:
             self.current_time = self.parameters["time"]
         else:
-            self.tau_frw, self.t_frw, self.dtau, self.n_frw, self.time_tot = friedman(
-                self.omega_matter,
-                self.omega_lambda,
-                1.0 - self.omega_matter - self.omega_lambda,
-            )
-
-            age = self.parameters["time"]
-            iage = 1 + int(10.0 * age / self.dtau)
-            iage = np.min([iage, self.n_frw // 2 + (iage - self.n_frw // 2) // 10])
-
-            try:
-                self.time_simu = self.t_frw[iage] * (age - self.tau_frw[iage - 1]) / (
-                    self.tau_frw[iage] - self.tau_frw[iage - 1]
-                ) + self.t_frw[iage - 1] * (age - self.tau_frw[iage]) / (
-                    self.tau_frw[iage - 1] - self.tau_frw[iage]
-                )
-
-                self.current_time = (
-                    (self.time_tot + self.time_simu)
-                    / (self.hubble_constant * 1e7 / cm_per_mpc)
-                    / self.parameters["unit_t"]
-                )
-            except IndexError:
-                mylog.warning(
-                    "Yt could not convert conformal time to physical time. "
-                    "Yt will assume the simulation is *not* cosmological."
-                )
-                self.cosmological_simulation = 0
-                self.current_time = self.parameters["time"]
+            aexp_grid = np.geomspace(1e-3, 1, 2_000, endpoint=False)
+            z_grid = 1 / aexp_grid - 1
+            self.tau_frw = tau_frw(self, z_grid)
+            self.t_frw = t_frw(self, z_grid)
+            self.current_time = t_frw(self, self.current_redshift).to("Gyr")
 
         if self.num_groups > 0:
             self.group_size = rheader["ncpu"] // self.num_groups
diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py
index 989d112860a..fcb88a0d441 100644
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -10,6 +10,7 @@
 from yt.fields.field_info_container import FieldInfoContainer
 from yt.frontends.ramses.io import convert_ramses_conformal_time_to_physical_age
 from yt.utilities.cython_fortran_utils import FortranFile
+from yt.utilities.lib.cosmology_time import t_frw
 from yt.utilities.linear_interpolators import BilinearFieldInterpolator
 from yt.utilities.logger import ytLogger as mylog
 from yt.utilities.physical_constants import (
@@ -179,21 +180,39 @@ class RAMSESFieldInfo(FieldInfoContainer):
     def setup_particle_fields(self, ptype):
         super().setup_particle_fields(ptype)
 
+        def star_age_from_conformal_cosmo(field, data):
+            conformal_age = data[ptype, "conformal_birth_time"]
+            birth_time = convert_ramses_conformal_time_to_physical_age(
+                data.ds, conformal_age
+            )
+            return data.ds.current_time - birth_time
+
+        def star_age_from_physical_cosmo(field, data):
+            H0 = float(
+                data.ds.quan(data.ds.hubble_constant * 100, "km/s/Mpc").to("1/Gyr")
+            )
+            times = data[ptype, "particle_birth_time"].value
+            time_tot = t_frw(data.ds, 0) * H0
+            birth_time = ((time_tot + times) / H0,)
+            return data.ds.apply_units(
+                data.ds.current_time.to("Gyr") - birth_time, "Gyr"
+            )
+
         def star_age(field, data):
-            if data.ds.cosmological_simulation and data.ds.use_conformal_time:
-                conformal_age = data[ptype, "conformal_birth_time"]
-                physical_age = convert_ramses_conformal_time_to_physical_age(
-                    data.ds, conformal_age
-                )
-                return data.ds.arr(physical_age, "code_time")
-            else:
-                formation_time = data[ptype, "particle_birth_time"]
-                return data.ds.current_time - formation_time
+            formation_time = data[ptype, "particle_birth_time"]
+            return data.ds.current_time - formation_time
+
+        if self.ds.cosmological_simulation and self.ds.use_conformal_time:
+            fun = star_age_from_conformal_cosmo
+        elif self.ds.cosmological_simulation:
+            fun = star_age_from_physical_cosmo
+        else:
+            fun = star_age
 
         self.add_field(
             (ptype, "star_age"),
             sampling_type="particle",
-            function=star_age,
+            function=fun,
             units=self.ds.unit_system["time"],
         )
 
diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index 0afa1464a9e..8e743c749c8 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -50,6 +50,12 @@ def convert_ramses_conformal_time_to_physical_age(
     physical_age : np.ndarray
         The physical age in code units
     """
+    h0 = ds.hubble_constant
+    tau_bins = ds.tau_frw * h0
+    t_bins = ds.t_frw
+
+    return ds.arr(np.interp(conformal_time, tau_bins, t_bins, right=0.0), t_bins.units)
+
     tf = ds.t_frw
     dtau = ds.dtau
     tauf = ds.tau_frw
diff --git a/yt/frontends/ramses/particle_handlers.py b/yt/frontends/ramses/particle_handlers.py
index 32da4d1aca3..dd9b34e173b 100644
--- a/yt/frontends/ramses/particle_handlers.py
+++ b/yt/frontends/ramses/particle_handlers.py
@@ -315,13 +315,13 @@ def build_iterator():
             )
             self.ds._warned_extra_fields["io"] = True
 
-        if self.has_birth_file:
-            # Sentinel value: this will prevent birth times from being read from
-            # the partXXXXX.outYYYYY files
-            field_offsets[ptype, "particle_birth_time"] = -1
-            _pfields[ptype, "particle_birth_time"] = "d"
-            self.ds.use_conformal_time = False
-        elif (
+        # if self.has_birth_file:
+        #     # Sentinel value: this will prevent birth times from being read from
+        #     # the partXXXXX.outYYYYY files
+        #     field_offsets[ptype, "particle_birth_time"] = -1
+        #     _pfields[ptype, "particle_birth_time"] = "d"
+        #     self.ds.use_conformal_time = False
+        if (
             self.ds.use_conformal_time
             and (ptype, "particle_birth_time") in field_offsets
         ):
diff --git a/yt/utilities/lib/cosmology_time.pyx b/yt/utilities/lib/cosmology_time.pyx
index 2a9d23cc723..90dbe87b71b 100644
--- a/yt/utilities/lib/cosmology_time.pyx
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -1,23 +1,92 @@
+# distutils: language = c++
 # distutils: libraries = STD_LIBS
+
 cimport numpy as np
 
 import numpy as np
 
 from libc.math cimport sqrt
+import cython
+from scipy.integrate import quad
+
+
+
+@cython.cdivision(True)
+cdef inline double _a_dot(double a, double h0, double om_m, double om_l) noexcept:
+    om_k = 1.0 - om_m - om_l
+    return h0 * a * sqrt(om_m * (a ** -3) + om_k * (a ** -2) + om_l)
+
+
+@cython.cdivision(True)
+cpdef double _a_dot_recip(double a, double h0, double om_m, double om_l):
+    return 1. / _a_dot(a, h0, om_m, om_l)
 
 
-cdef double dadtau(double aexp_tau,double O_mat_0,double O_vac_0,double O_k_0):
+cdef inline double _da_dtau(double a, double h0, double om_m, double om_l) noexcept:
+    return a**2 * _a_dot(a, h0, om_m, om_l)
+
+@cython.cdivision(True)
+cpdef double _da_dtau_recip(double a, double h0, double om_m, double om_l) noexcept:
+    return 1. / _da_dtau(a, h0, om_m, om_l)
+
+
+cdef double dadtau(double aexp_tau,double O_mat_0,double O_vac_0,double O_k_0) noexcept:
     cdef double aexp_tau3 = aexp_tau * aexp_tau * aexp_tau
     return sqrt( aexp_tau3 * (O_mat_0 + O_vac_0*aexp_tau3 + O_k_0*aexp_tau) )
 
-cdef double dadt(double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
+@cython.cdivision(True)
+cdef double dadt(double aexp_t,double O_mat_0,double O_vac_0,double O_k_0) noexcept:
     cdef double aexp_t3 = aexp_t * aexp_t * aexp_t
     return sqrt( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t3 + O_k_0*aexp_t) )
 
 
-cdef step_cosmo(double alpha,double tau,double aexp_tau,double t,double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
-    cdef double dtau,aexp_tau_pre,dt,aexp_t_pre
-
+def t_frw(ds, z):
+    from scipy.integrate import quad
+    aexp = 1 / (1 + z)
+
+    h0 = ds.hubble_constant
+    om_m = ds.omega_matter
+    om_l = ds.omega_lambda
+    conv = ds.quan(0.01, "Mpc/km*s").to("Gyr")
+
+    if isinstance(z, float):
+        return ds.quan(
+            quad(_a_dot_recip, 0, aexp, args=(h0, om_m, om_l))[0],
+            units=conv,
+        )
+
+
+    return ds.arr(
+        [quad(_a_dot_recip, 0, a, args=(h0, om_m, om_l))[0] for a in aexp],
+        units=conv,
+    )
+
+def tau_frw(ds, z):
+    from scipy.integrate import quad
+    aexp = 1 / (1 + z)
+
+    h0 = ds.hubble_constant
+    om_m = ds.omega_matter
+    om_l = ds.omega_lambda
+
+    if isinstance(z, float):
+        return quad(_da_dtau_recip, 1, aexp, args=(h0, om_m, om_l))[0]
+
+    return np.asarray(
+        [quad(_da_dtau_recip, 1, a, args=(h0, om_m, om_l))[0] for a in aexp],
+    )
+
+@cython.cdivision(True)
+cdef void step_cosmo(
+    const double alpha,
+    double &tau,
+    double &aexp_tau,
+    double &t,
+    double &aexp_t,
+    const double O_mat_0,
+    const double O_vac_0,
+    const double O_k_0,
+) noexcept:
     dtau = alpha * aexp_tau / dadtau(aexp_tau,O_mat_0,O_vac_0,O_k_0)
     aexp_tau_pre = aexp_tau - dadtau(aexp_tau,O_mat_0,O_vac_0,O_k_0)*dtau/2.0
     aexp_tau = aexp_tau - dadtau(aexp_tau_pre,O_mat_0,O_vac_0,O_k_0)*dtau
@@ -28,18 +97,26 @@ cdef step_cosmo(double alpha,double tau,double aexp_tau,double t,double aexp_t,d
     aexp_t = aexp_t - dadt(aexp_t_pre,O_mat_0,O_vac_0,O_k_0)*dt
     t = t - dt
 
-    return tau,aexp_tau,t,aexp_t
-
-
 cpdef friedman(double O_mat_0,double O_vac_0,double O_k_0):
-    cdef double alpha=1.e-5,aexp_min=1.e-3,aexp_tau=1.,aexp_t=1.,tau=0.,t=0.
-    cdef int nstep=0,ntable=1000,n_out
-    cdef np.ndarray[double,mode='c'] t_out=np.zeros([ntable+1]),tau_out=np.zeros([ntable+1])
-    cdef double age_tot,delta_tau,next_tau
+    cdef double alpha=1e-4 ,aexp_min=1e-3, aexp_tau=1, aexp_t=1, tau=0, t=0
+    cdef int nstep=0, ntable=1000, n_out
+    cdef np.ndarray[double,mode='c'] t_out = np.zeros([ntable+1])
+    cdef np.ndarray[double,mode='c'] tau_out = np.zeros([ntable+1])
+    cdef double age_tot, delta_tau, next_tau
 
     while aexp_tau >= aexp_min or aexp_t >= aexp_min:
-       nstep = nstep + 1
-       tau,aexp_tau,t,aexp_t = step_cosmo(alpha,tau,aexp_tau,t,aexp_t,O_mat_0,O_vac_0,O_k_0)
+        nstep = nstep + 1
+
+        step_cosmo(
+            alpha,
+            tau,
+            aexp_tau,
+            t,
+            aexp_t,
+            O_mat_0,
+            O_vac_0,
+            O_k_0,
+        )
 
     age_tot=-t
     if nstep < ntable :
@@ -60,7 +137,16 @@ cpdef friedman(double O_mat_0,double O_vac_0,double O_k_0):
     next_tau = tau + delta_tau/10.
 
     while n_out < ntable/2 :
-        tau,aexp_tau,t,aexp_t = step_cosmo(alpha,tau,aexp_tau,t,aexp_t,O_mat_0,O_vac_0,O_k_0)
+        step_cosmo(
+            alpha,
+            tau,
+            aexp_tau,
+            t,
+            aexp_t,
+            O_mat_0,
+            O_vac_0,
+            O_k_0,
+        )
 
         if tau < next_tau:
             n_out = n_out + 1
@@ -69,7 +155,16 @@ cpdef friedman(double O_mat_0,double O_vac_0,double O_k_0):
             next_tau = next_tau + delta_tau/10.
 
     while aexp_tau >= aexp_min or aexp_t >= aexp_min:
-        tau,aexp_tau,t,aexp_t = step_cosmo(alpha,tau,aexp_tau,t,aexp_t,O_mat_0,O_vac_0,O_k_0)
+        step_cosmo(
+            alpha,
+            tau,
+            aexp_tau,
+            t,
+            aexp_t,
+            O_mat_0,
+            O_vac_0,
+            O_k_0,
+        )
 
         if tau < next_tau:
             n_out = n_out + 1

From 01098121eb9a0dd71c7baa20aae454f3e43d4469 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@cphyc.me>
Date: Thu, 2 Nov 2023 17:28:23 +0000
Subject: [PATCH 03/13] Also work for non-conformal simulations

---
 yt/frontends/ramses/fields.py       | 11 +++++------
 yt/utilities/lib/cosmology_time.pyx |  5 ++---
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py
index fcb88a0d441..3624224ba90 100644
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -191,12 +191,11 @@ def star_age_from_physical_cosmo(field, data):
             H0 = float(
                 data.ds.quan(data.ds.hubble_constant * 100, "km/s/Mpc").to("1/Gyr")
             )
-            times = data[ptype, "particle_birth_time"].value
-            time_tot = t_frw(data.ds, 0) * H0
-            birth_time = ((time_tot + times) / H0,)
-            return data.ds.apply_units(
-                data.ds.current_time.to("Gyr") - birth_time, "Gyr"
-            )
+            times = data[ptype, "conformal_birth_time"].value
+            time_tot = float(t_frw(data.ds, 0) * H0)
+            birth_time = (time_tot + times) / H0
+            t_out = float(data.ds.current_time.to("Gyr"))
+            return data.apply_units(t_out - birth_time, "Gyr")
 
         def star_age(field, data):
             formation_time = data[ptype, "particle_birth_time"]
diff --git a/yt/utilities/lib/cosmology_time.pyx b/yt/utilities/lib/cosmology_time.pyx
index 90dbe87b71b..0a2eb9d0ec7 100644
--- a/yt/utilities/lib/cosmology_time.pyx
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -49,13 +49,12 @@ def t_frw(ds, z):
     om_l = ds.omega_lambda
     conv = ds.quan(0.01, "Mpc/km*s").to("Gyr")
 
-    if isinstance(z, float):
+    if isinstance(z, (int, float)):
         return ds.quan(
             quad(_a_dot_recip, 0, aexp, args=(h0, om_m, om_l))[0],
             units=conv,
         )
 
-
     return ds.arr(
         [quad(_a_dot_recip, 0, a, args=(h0, om_m, om_l))[0] for a in aexp],
         units=conv,
@@ -69,7 +68,7 @@ def tau_frw(ds, z):
     om_m = ds.omega_matter
     om_l = ds.omega_lambda
 
-    if isinstance(z, float):
+    if isinstance(z, (int, float)):
         return quad(_da_dtau_recip, 1, aexp, args=(h0, om_m, om_l))[0]
 
     return np.asarray(

From 1a7123332f2f48719dcad1b2c2508e78a0f41212 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@cphyc.me>
Date: Thu, 2 Nov 2023 17:30:44 +0000
Subject: [PATCH 04/13] Cleanup cosmology_time

---
 yt/utilities/lib/cosmology_time.pyx | 111 ----------------------------
 1 file changed, 111 deletions(-)

diff --git a/yt/utilities/lib/cosmology_time.pyx b/yt/utilities/lib/cosmology_time.pyx
index 0a2eb9d0ec7..8e2ba37fdc6 100644
--- a/yt/utilities/lib/cosmology_time.pyx
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -30,15 +30,6 @@ cpdef double _da_dtau_recip(double a, double h0, double om_m, double om_l) noexc
     return 1. / _da_dtau(a, h0, om_m, om_l)
 
 
-cdef double dadtau(double aexp_tau,double O_mat_0,double O_vac_0,double O_k_0) noexcept:
-    cdef double aexp_tau3 = aexp_tau * aexp_tau * aexp_tau
-    return sqrt( aexp_tau3 * (O_mat_0 + O_vac_0*aexp_tau3 + O_k_0*aexp_tau) )
-
-@cython.cdivision(True)
-cdef double dadt(double aexp_t,double O_mat_0,double O_vac_0,double O_k_0) noexcept:
-    cdef double aexp_t3 = aexp_t * aexp_t * aexp_t
-    return sqrt( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t3 + O_k_0*aexp_t) )
-
 
 def t_frw(ds, z):
     from scipy.integrate import quad
@@ -74,105 +65,3 @@ def tau_frw(ds, z):
     return np.asarray(
         [quad(_da_dtau_recip, 1, a, args=(h0, om_m, om_l))[0] for a in aexp],
     )
-
-@cython.cdivision(True)
-cdef void step_cosmo(
-    const double alpha,
-    double &tau,
-    double &aexp_tau,
-    double &t,
-    double &aexp_t,
-    const double O_mat_0,
-    const double O_vac_0,
-    const double O_k_0,
-) noexcept:
-    dtau = alpha * aexp_tau / dadtau(aexp_tau,O_mat_0,O_vac_0,O_k_0)
-    aexp_tau_pre = aexp_tau - dadtau(aexp_tau,O_mat_0,O_vac_0,O_k_0)*dtau/2.0
-    aexp_tau = aexp_tau - dadtau(aexp_tau_pre,O_mat_0,O_vac_0,O_k_0)*dtau
-    tau = tau - dtau
-
-    dt = alpha * aexp_t / dadt(aexp_t,O_mat_0,O_vac_0,O_k_0)
-    aexp_t_pre = aexp_t - dadt(aexp_t,O_mat_0,O_vac_0,O_k_0)*dt/2.0
-    aexp_t = aexp_t - dadt(aexp_t_pre,O_mat_0,O_vac_0,O_k_0)*dt
-    t = t - dt
-
-cpdef friedman(double O_mat_0,double O_vac_0,double O_k_0):
-    cdef double alpha=1e-4 ,aexp_min=1e-3, aexp_tau=1, aexp_t=1, tau=0, t=0
-    cdef int nstep=0, ntable=1000, n_out
-    cdef np.ndarray[double,mode='c'] t_out = np.zeros([ntable+1])
-    cdef np.ndarray[double,mode='c'] tau_out = np.zeros([ntable+1])
-    cdef double age_tot, delta_tau, next_tau
-
-    while aexp_tau >= aexp_min or aexp_t >= aexp_min:
-        nstep = nstep + 1
-
-        step_cosmo(
-            alpha,
-            tau,
-            aexp_tau,
-            t,
-            aexp_t,
-            O_mat_0,
-            O_vac_0,
-            O_k_0,
-        )
-
-    age_tot=-t
-    if nstep < ntable :
-        ntable = nstep
-        alpha = alpha / 2.
-
-    delta_tau = 20.*tau/ntable/11.
-
-    aexp_tau = 1.
-    aexp_t = 1.
-    tau = 0.
-    t = 0.
-
-    n_out = 0
-    t_out[n_out] = t
-    tau_out[n_out] = tau
-
-    next_tau = tau + delta_tau/10.
-
-    while n_out < ntable/2 :
-        step_cosmo(
-            alpha,
-            tau,
-            aexp_tau,
-            t,
-            aexp_t,
-            O_mat_0,
-            O_vac_0,
-            O_k_0,
-        )
-
-        if tau < next_tau:
-            n_out = n_out + 1
-            t_out[n_out] = t
-            tau_out[n_out] = tau
-            next_tau = next_tau + delta_tau/10.
-
-    while aexp_tau >= aexp_min or aexp_t >= aexp_min:
-        step_cosmo(
-            alpha,
-            tau,
-            aexp_tau,
-            t,
-            aexp_t,
-            O_mat_0,
-            O_vac_0,
-            O_k_0,
-        )
-
-        if tau < next_tau:
-            n_out = n_out + 1
-            t_out[n_out] = t
-            tau_out[n_out] = tau
-            next_tau = next_tau + delta_tau
-
-    n_out = ntable
-    t_out[n_out] = t
-    tau_out[n_out] = tau
-
-    return tau_out,t_out,delta_tau,ntable,age_tot

From 68cb79e82f5368f58e412cc90ffe40019be76237 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@fysik.lu.se>
Date: Fri, 17 May 2024 16:59:49 +0200
Subject: [PATCH 05/13] Convert to code time

---
 yt/frontends/ramses/io.py                | 3 ++-
 yt/frontends/ramses/particle_handlers.py | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index 8e743c749c8..b90e2dbc810 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -15,6 +15,7 @@
 from yt.utilities.io_handler import BaseIOHandler
 from yt.utilities.logger import ytLogger as mylog
 from yt.utilities.physical_ratios import cm_per_km, cm_per_mpc
+from unyt import unyt_array
 
 if TYPE_CHECKING:
     import os
@@ -34,7 +35,7 @@ def convert_ramses_ages(ds, conformal_ages):
 
 def convert_ramses_conformal_time_to_physical_age(
     ds, conformal_time: np.ndarray
-) -> np.ndarray:
+) -> unyt_array:
     """
     Convert conformal times (as defined in RAMSES) to physical age.
 
diff --git a/yt/frontends/ramses/particle_handlers.py b/yt/frontends/ramses/particle_handlers.py
index dd9b34e173b..8303fe2b0a1 100644
--- a/yt/frontends/ramses/particle_handlers.py
+++ b/yt/frontends/ramses/particle_handlers.py
@@ -370,7 +370,7 @@ def handle_field(
         # Otherwise, convert conformal time to physical age
         ds = self.ds
         conformal_time = data_dict[field]
-        physical_age = convert_ramses_conformal_time_to_physical_age(ds, conformal_time)
+        physical_age = convert_ramses_conformal_time_to_physical_age(ds, conformal_time).to("code_time").v
         current_time = ds.current_time.in_units("code_time").v
         # arbitrarily set particles with zero conformal_age to zero
         # particle_age. This corresponds to DM particles.

From 68fd745967ab97367c2ab8c1dcdcd476b24b4f70 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@fysik.lu.se>
Date: Fri, 17 May 2024 17:04:21 +0200
Subject: [PATCH 06/13] Interp to min/max time when running out of bounds

---
 yt/frontends/ramses/io.py | 32 ++++++--------------------------
 1 file changed, 6 insertions(+), 26 deletions(-)

diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index b90e2dbc810..3f37976a0c3 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -3,6 +3,7 @@
 from typing import TYPE_CHECKING, Union
 
 import numpy as np
+from unyt import unyt_array
 
 from yt._maintenance.deprecation import issue_deprecation_warning
 from yt.frontends.ramses.definitions import VAR_DESC_RE, VERSION_RE
@@ -14,8 +15,6 @@
 )
 from yt.utilities.io_handler import BaseIOHandler
 from yt.utilities.logger import ytLogger as mylog
-from yt.utilities.physical_ratios import cm_per_km, cm_per_mpc
-from unyt import unyt_array
 
 if TYPE_CHECKING:
     import os
@@ -55,31 +54,12 @@ def convert_ramses_conformal_time_to_physical_age(
     tau_bins = ds.tau_frw * h0
     t_bins = ds.t_frw
 
-    return ds.arr(np.interp(conformal_time, tau_bins, t_bins, right=0.0), t_bins.units)
-
-    tf = ds.t_frw
-    dtau = ds.dtau
-    tauf = ds.tau_frw
-    h100 = ds.hubble_constant
-    nOver2 = ds.n_frw / 2
-    unit_t = ds.parameters["unit_t"]
-    tsim = ds.time_simu
-
-    t_scale = 1.0 / (h100 * 100 * cm_per_km / cm_per_mpc) / unit_t
-
-    # calculate index into lookup table (n_frw elements in
-    # lookup table)
-    dage = 1 + (10 * conformal_time / dtau)
-    dage = np.minimum(dage, nOver2 + (dage - nOver2) / 10.0)
-    iage = np.array(dage, dtype=np.int32)
-
-    # linearly interpolate physical times from tf and tauf lookup
-    # tables.
-    t = tf[iage] * (conformal_time - tauf[iage - 1]) / (tauf[iage] - tauf[iage - 1])
-    t = t + (
-        tf[iage - 1] * (conformal_time - tauf[iage]) / (tauf[iage - 1] - tauf[iage])
+    return ds.arr(
+        np.interp(
+            conformal_time, tau_bins, t_bins, right=t_bins.max(), left=t_bins.min()
+        ),
+        t_bins.units,
     )
-    return (tsim - t) * t_scale
 
 
 def _ramses_particle_binary_file_handler(particle_handler, subset, fields, count):

From 74940505966c5f1ea7ab0bde7484857c6ffd66ae Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@fysik.lu.se>
Date: Fri, 17 May 2024 17:19:31 +0200
Subject: [PATCH 07/13] Fix linting

---
 yt/frontends/ramses/particle_handlers.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/yt/frontends/ramses/particle_handlers.py b/yt/frontends/ramses/particle_handlers.py
index 8303fe2b0a1..a295d66217d 100644
--- a/yt/frontends/ramses/particle_handlers.py
+++ b/yt/frontends/ramses/particle_handlers.py
@@ -370,7 +370,11 @@ def handle_field(
         # Otherwise, convert conformal time to physical age
         ds = self.ds
         conformal_time = data_dict[field]
-        physical_age = convert_ramses_conformal_time_to_physical_age(ds, conformal_time).to("code_time").v
+        physical_age = (
+            convert_ramses_conformal_time_to_physical_age(ds, conformal_time)
+            .to("code_time")
+            .v
+        )
         current_time = ds.current_time.in_units("code_time").v
         # arbitrarily set particles with zero conformal_age to zero
         # particle_age. This corresponds to DM particles.

From 1e6055e75e578da0ed14283bd8f2dbe16a818e30 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@fysik.lu.se>
Date: Fri, 17 May 2024 17:30:48 +0200
Subject: [PATCH 08/13] Require scipy for RAMSES only

---
 pyproject.toml                      | 2 +-
 yt/utilities/lib/cosmology_time.pyx | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 98656c17d0f..c7e7d1b82dc 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -121,7 +121,7 @@ open-pmd = ["yt[HDF5]"]
 owls = ["yt[HDF5]"]
 owls-subfind = ["yt[HDF5]"]
 parthenon = ["yt[HDF5]"]
-ramses = ["yt[Fortran]"]
+ramses = ["yt[Fortran]", "scipy"]
 rockstar = []
 sdf = ["requests>=2.20.0"]
 stream = []
diff --git a/yt/utilities/lib/cosmology_time.pyx b/yt/utilities/lib/cosmology_time.pyx
index 8e2ba37fdc6..9d1d65c3fb0 100644
--- a/yt/utilities/lib/cosmology_time.pyx
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -7,7 +7,6 @@ import numpy as np
 
 from libc.math cimport sqrt
 import cython
-from scipy.integrate import quad
 
 
 

From 08ff40d67bf4b6c8de587018fe0882b4ff6e85b2 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@cphyc.me>
Date: Thu, 23 May 2024 16:57:52 +0200
Subject: [PATCH 09/13] Clip min/max times to 0, tsim

---
 yt/frontends/ramses/io.py                 |  9 +++++++--
 yt/frontends/ramses/tests/test_outputs.py | 10 ++++++----
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index 3f37976a0c3..66e72130023 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -54,9 +54,14 @@ def convert_ramses_conformal_time_to_physical_age(
     tau_bins = ds.tau_frw * h0
     t_bins = ds.t_frw
 
+    min_time = 0
+    max_time = ds.current_time.to(t_bins.units)
+
     return ds.arr(
-        np.interp(
-            conformal_time, tau_bins, t_bins, right=t_bins.max(), left=t_bins.min()
+        np.clip(
+            np.interp(conformal_time, tau_bins, t_bins, right=max_time, left=min_time),
+            min_time,
+            max_time.value,
         ),
         t_bins.units,
     )
diff --git a/yt/frontends/ramses/tests/test_outputs.py b/yt/frontends/ramses/tests/test_outputs.py
index 064c57b78f8..bc76f27d7e6 100644
--- a/yt/frontends/ramses/tests/test_outputs.py
+++ b/yt/frontends/ramses/tests/test_outputs.py
@@ -119,11 +119,13 @@ def test_unit_cosmo():
         # expected_raw_time = 1.119216564055017 # in ramses unit
         # expected_time = 3.756241729312462e17 # in seconds
 
-        expected_raw_time = 1.121279694787743  # in ramses unit
-        assert_equal(ds.current_time.value, expected_raw_time)
+        expected_raw_time = 1.1213297131656712  # in ramses unit
+        np.testing.assert_equal(
+            ds.current_time.to("code_time").value, expected_raw_time
+        )
 
-        expected_time = 3.7631658742904595e17  # in seconds
-        assert_equal(ds.current_time.in_units("s").value, expected_time)
+        expected_time = 3.7633337427123904e17  # in seconds
+        np.testing.assert_equal(ds.current_time.to("s").value, expected_time)
 
 
 ramsesExtraFieldsSmall = "ramses_extra_fields_small/output_00001"

From c37f115b9deda32344e3faf00c3e03ff533e339c Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@cphyc.me>
Date: Fri, 24 May 2024 11:23:52 +0200
Subject: [PATCH 10/13] Strip unit before clipping/interpolating

---
 yt/frontends/ramses/io.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index 66e72130023..38312e35f4a 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -59,7 +59,13 @@ def convert_ramses_conformal_time_to_physical_age(
 
     return ds.arr(
         np.clip(
-            np.interp(conformal_time, tau_bins, t_bins, right=max_time, left=min_time),
+            np.interp(
+                conformal_time.value,
+                tau_bins,
+                t_bins.value,
+                right=max_time,
+                left=min_time,
+            ),
             min_time,
             max_time.value,
         ),

From 9b6fbf6a81ae3bb087b830bf2f55d3c99b92f67c Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <corentin.cadiou@cphyc.me>
Date: Fri, 24 May 2024 11:58:50 +0200
Subject: [PATCH 11/13] Rename functions for better good

---
 yt/frontends/ramses/fields.py            |  4 ++--
 yt/frontends/ramses/io.py                |  8 ++++----
 yt/frontends/ramses/particle_handlers.py | 16 ++++++++--------
 3 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/yt/frontends/ramses/fields.py b/yt/frontends/ramses/fields.py
index 3624224ba90..4f8ddea9261 100644
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -8,7 +8,7 @@
 from yt._typing import KnownFieldsT
 from yt.fields.field_detector import FieldDetector
 from yt.fields.field_info_container import FieldInfoContainer
-from yt.frontends.ramses.io import convert_ramses_conformal_time_to_physical_age
+from yt.frontends.ramses.io import convert_ramses_conformal_time_to_physical_time
 from yt.utilities.cython_fortran_utils import FortranFile
 from yt.utilities.lib.cosmology_time import t_frw
 from yt.utilities.linear_interpolators import BilinearFieldInterpolator
@@ -182,7 +182,7 @@ def setup_particle_fields(self, ptype):
 
         def star_age_from_conformal_cosmo(field, data):
             conformal_age = data[ptype, "conformal_birth_time"]
-            birth_time = convert_ramses_conformal_time_to_physical_age(
+            birth_time = convert_ramses_conformal_time_to_physical_time(
                 data.ds, conformal_age
             )
             return data.ds.current_time - birth_time
diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index 38312e35f4a..c4f0960e032 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -29,14 +29,14 @@ def convert_ramses_ages(ds, conformal_ages):
         stacklevel=3,
         since="4.0.3",
     )
-    return convert_ramses_conformal_time_to_physical_age(ds, conformal_ages)
+    return convert_ramses_conformal_time_to_physical_time(ds, conformal_ages)
 
 
-def convert_ramses_conformal_time_to_physical_age(
+def convert_ramses_conformal_time_to_physical_time(
     ds, conformal_time: np.ndarray
 ) -> unyt_array:
     """
-    Convert conformal times (as defined in RAMSES) to physical age.
+    Convert conformal times (as defined in RAMSES) to physical times.
 
     Arguments
     ---------
@@ -60,7 +60,7 @@ def convert_ramses_conformal_time_to_physical_age(
     return ds.arr(
         np.clip(
             np.interp(
-                conformal_time.value,
+                conformal_time,
                 tau_bins,
                 t_bins.value,
                 right=max_time,
diff --git a/yt/frontends/ramses/particle_handlers.py b/yt/frontends/ramses/particle_handlers.py
index a295d66217d..d151dea722c 100644
--- a/yt/frontends/ramses/particle_handlers.py
+++ b/yt/frontends/ramses/particle_handlers.py
@@ -19,7 +19,7 @@
     _ramses_particle_csv_file_handler,
     _read_part_binary_file_descriptor,
     _read_part_csv_file_descriptor,
-    convert_ramses_conformal_time_to_physical_age,
+    convert_ramses_conformal_time_to_physical_time,
 )
 
 if TYPE_CHECKING:
@@ -370,15 +370,14 @@ def handle_field(
         # Otherwise, convert conformal time to physical age
         ds = self.ds
         conformal_time = data_dict[field]
-        physical_age = (
-            convert_ramses_conformal_time_to_physical_age(ds, conformal_time)
+        physical_time = (
+            convert_ramses_conformal_time_to_physical_time(ds, conformal_time)
             .to("code_time")
             .v
         )
-        current_time = ds.current_time.in_units("code_time").v
         # arbitrarily set particles with zero conformal_age to zero
         # particle_age. This corresponds to DM particles.
-        data_dict[field] = np.where(conformal_time > 0, current_time - physical_age, 0)
+        data_dict[field] = np.where(conformal_time != 0, physical_time, 0)
 
 
 class SinkParticleFileHandler(ParticleFileHandler):
@@ -509,8 +508,9 @@ def handle_field(
         # onvert conformal time to physical age
         ds = self.ds
         conformal_time = data_dict[field]
-        physical_age = convert_ramses_conformal_time_to_physical_age(ds, conformal_time)
-        current_time = ds.current_time.in_units("code_time").v
+        physical_time = convert_ramses_conformal_time_to_physical_time(
+            ds, conformal_time
+        )
         # arbitrarily set particles with zero conformal_age to zero
         # particle_age. This corresponds to DM particles.
-        data_dict[field] = np.where(conformal_time > 0, current_time - physical_age, 0)
+        data_dict[field] = np.where(conformal_time > 0, physical_time, 0)

From 2e7641228998f34c63bf9fb801d953d7f4fec7d8 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <contact@cphyc.me>
Date: Thu, 3 Oct 2024 13:10:57 +0200
Subject: [PATCH 12/13] Apply suggestions from code review

Co-authored-by: Chris Havlin <chris.havlin@gmail.com>
---
 yt/frontends/ramses/io.py                | 2 +-
 yt/frontends/ramses/particle_handlers.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/yt/frontends/ramses/io.py b/yt/frontends/ramses/io.py
index c4f0960e032..6f241631a8a 100644
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -110,7 +110,7 @@ def _ramses_particle_binary_file_handler(particle_handler, subset, fields, count
             if field[1].startswith("particle_position"):
                 np.divide(tr[field], ds["boxlen"], tr[field])
 
-            # Handle over to field handler for special cases, like particle_birth_times
+            # Hand over to field handler for special cases, like particle_birth_times
             particle_handler.handle_field(field, tr)
 
     return tr
diff --git a/yt/frontends/ramses/particle_handlers.py b/yt/frontends/ramses/particle_handlers.py
index d151dea722c..3cf181cb1fd 100644
--- a/yt/frontends/ramses/particle_handlers.py
+++ b/yt/frontends/ramses/particle_handlers.py
@@ -505,7 +505,7 @@ def handle_field(
         if not (fname == "particle_birth_time" and self.ds.cosmological_simulation):
             return
 
-        # onvert conformal time to physical age
+        # convert conformal time to physical age
         ds = self.ds
         conformal_time = data_dict[field]
         physical_time = convert_ramses_conformal_time_to_physical_time(

From 2769152ece5d745187912c4573279897df09ec01 Mon Sep 17 00:00:00 2001
From: Corentin Cadiou <contact@cphyc.me>
Date: Thu, 3 Oct 2024 13:11:46 +0200
Subject: [PATCH 13/13] Remove commented-out lines

---
 yt/frontends/ramses/particle_handlers.py | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/yt/frontends/ramses/particle_handlers.py b/yt/frontends/ramses/particle_handlers.py
index 3cf181cb1fd..ff963e578bf 100644
--- a/yt/frontends/ramses/particle_handlers.py
+++ b/yt/frontends/ramses/particle_handlers.py
@@ -315,12 +315,6 @@ def build_iterator():
             )
             self.ds._warned_extra_fields["io"] = True
 
-        # if self.has_birth_file:
-        #     # Sentinel value: this will prevent birth times from being read from
-        #     # the partXXXXX.outYYYYY files
-        #     field_offsets[ptype, "particle_birth_time"] = -1
-        #     _pfields[ptype, "particle_birth_time"] = "d"
-        #     self.ds.use_conformal_time = False
         if (
             self.ds.use_conformal_time
             and (ptype, "particle_birth_time") in field_offsets