diff --git a/src/osipi/_aif.py b/src/osipi/_aif.py
index 1af666f..bbb995b 100644
--- a/src/osipi/_aif.py
+++ b/src/osipi/_aif.py
@@ -1,14 +1,14 @@
 import numpy as np
-from numpy.typing import NDArray
+from numpy.typing import ArrayLike, NDArray
 
 
 def aif_parker(
-    t: NDArray[np.floating], BAT: np.floating = 0.0, Hct: np.floating = 0.0
+    t: ArrayLike, BAT: np.floating = 0.0, Hct: np.floating = 0.0
 ) -> NDArray[np.floating]:
     """AIF model as defined by Parker et al (2005)
 
     Args:
-        t (NDArray[np.floating]): array of time points in units of sec. [OSIPI code Q.GE1.004]
+        t (ArrayLike): array of time points in units of sec. [OSIPI code Q.GE1.004]
         BAT (np.floating, optional):
             Time in seconds before the bolus arrives. Defaults to 0. [OSIPI code Q.BA1.001]
         Hct (np.floating, optional):
@@ -45,8 +45,14 @@ def aif_parker(
         >>> plt.show()
 
     """
+    # Convert input to numpy array with appropriate dtype
+    t_arr = np.asarray(t)
+    # Ensure floating-point dtype
+    if not np.issubdtype(t_arr.dtype, np.floating):
+        t_arr = t_arr.astype(np.float64)
+
     # Convert from OSIPI units (sec) to units used internally (mins)
-    t_min = t / 60
+    t_min = t_arr / 60
     bat_min = BAT / 60
 
     t_offset = t_min - bat_min
@@ -71,7 +77,7 @@ def aif_parker(
     return pop_aif
 
 
-def aif_georgiou(t: NDArray[np.floating], BAT: np.floating = 0.0) -> NDArray[np.floating]:
+def aif_georgiou(t: ArrayLike, BAT: np.floating = 0.0) -> NDArray[np.floating]:
     """AIF model as defined by Georgiou et al.
 
     Note:
@@ -80,7 +86,7 @@ def aif_georgiou(t: NDArray[np.floating], BAT: np.floating = 0.0) -> NDArray[np.
         so nobody ever has to write this function again!
 
     Args:
-        t (NDArray[np.floating]): array of time points in units of sec. [OSIPI code Q.GE1.004]
+        t (ArrayLike): array of time points in units of sec. [OSIPI code Q.GE1.004]
         BAT (np.floating, optional):
             Time in seconds before the bolus arrives. Defaults to 0sec. [OSIPI code Q.BA1.001]
 
@@ -119,13 +125,12 @@ def aif_georgiou(t: NDArray[np.floating], BAT: np.floating = 0.0) -> NDArray[np.
 
     msg = "This function is not yet implemented \n"
     msg += (
-        "If you implement it yourself, please consider submitting it"
-        " as an OSIPI code contribution"
+        "If you implement it yourself, please consider submitting it as an OSIPI code contribution"
     )
     raise NotImplementedError(msg)
 
 
-def aif_weinmann(t: NDArray[np.floating], BAT: np.floating = 0.0) -> NDArray[np.floating]:
+def aif_weinmann(t: ArrayLike, BAT: np.floating = 0.0) -> NDArray[np.floating]:
     """AIF model as defined by Weinmann et al.
 
     Note:
@@ -134,7 +139,7 @@ def aif_weinmann(t: NDArray[np.floating], BAT: np.floating = 0.0) -> NDArray[np.
         so nobody ever has to write this function again!
 
     Args:
-        t (NDArray[np.floating]): array of time points in units of sec. [OSIPI code Q.GE1.004]
+        t (ArrayLike): array of time points in units of sec. [OSIPI code Q.GE1.004]
         BAT (np.floating, optional):
             Time in seconds before the bolus arrives. Defaults to 0sec. [OSIPI code Q.BA1.001]
 
@@ -172,7 +177,6 @@ def aif_weinmann(t: NDArray[np.floating], BAT: np.floating = 0.0) -> NDArray[np.
 
     msg = "This function is not yet implemented \n"
     msg += (
-        "If you implement it yourself, please consider submitting it"
-        " as an OSIPI code contribution"
+        "If you implement it yourself, please consider submitting it as an OSIPI code contribution"
     )
     raise NotImplementedError(msg)
diff --git a/src/osipi/_convolution.py b/src/osipi/_convolution.py
index 6bcccea..063e7f3 100755
--- a/src/osipi/_convolution.py
+++ b/src/osipi/_convolution.py
@@ -1,34 +1,42 @@
 import numpy as np
-from numpy.typing import NDArray
+from numpy.typing import ArrayLike, NDArray
 
 
-def exp_conv(
-    T: np.floating, t: NDArray[np.floating], a: NDArray[np.floating]
-) -> NDArray[np.floating]:
+def exp_conv(T: np.floating, t: ArrayLike, a: ArrayLike) -> NDArray[np.floating]:
     """Exponential convolution operation of (1/T)exp(-t/T) with a.
 
     Args:
         T (np.floating): exponent in time units
-        t (NDArray[np.floating]): array of time points
-        a (NDArray[np.floating]): array to be convolved with time exponential
+        t (ArrayLike): array of time points
+        a (ArrayLike): array to be convolved with time exponential
 
     Returns:
         NDArray[np.floating]: convolved array
     """
+    # Convert inputs to numpy arrays
+    t_arr = np.asarray(t)
+    a_arr = np.asarray(a)
+
+    # Ensure floating-point dtype while preserving original precision
+    if not np.issubdtype(t_arr.dtype, np.floating):
+        t_arr = t_arr.astype(np.float64)
+    if not np.issubdtype(a_arr.dtype, np.floating):
+        a_arr = a_arr.astype(np.float64)
+
     if T == 0:
-        return a
+        return a_arr
 
-    n = len(t)
-    f = np.zeros((n,))
+    n = len(t_arr)
+    f = np.zeros((n,), dtype=t_arr.dtype)  # Use same dtype as t_arr
 
-    x = (t[1 : n - 1] - t[0 : n - 2]) / T
-    da = (a[1 : n - 1] - a[0 : n - 2]) / x
+    x = (t_arr[1 : n - 1] - t_arr[0 : n - 2]) / T
+    da = (a_arr[1 : n - 1] - a_arr[0 : n - 2]) / x
 
     E = np.exp(-x)
     E0 = 1 - E
     E1 = x - E0
 
-    add = a[0 : n - 2] * E0 + da * E1
+    add = a_arr[0 : n - 2] * E0 + da * E1
 
     for i in range(0, n - 2):
         f[i + 1] = E[i] * f[i] + add[i]
diff --git a/src/osipi/_electromagnetic_property.py b/src/osipi/_electromagnetic_property.py
index 693ffd8..604ab91 100644
--- a/src/osipi/_electromagnetic_property.py
+++ b/src/osipi/_electromagnetic_property.py
@@ -1,9 +1,9 @@
 import numpy as np
-from numpy.typing import NDArray
+from numpy.typing import ArrayLike, NDArray
 
 
 def R1_to_C_linear_relaxivity(
-    R1: NDArray[np.floating], R10: np.floating, r1: np.floating
+    R1: ArrayLike, R10: np.floating, r1: np.floating
 ) -> NDArray[np.floating]:
     """
     Electromagnetic property inverse model:
@@ -12,7 +12,7 @@ def R1_to_C_linear_relaxivity(
     Converts R1 to tissue concentration
 
     Args:
-        R1 (1D array of np.floating):
+        R1 (ArrayLike):
             Vector of longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.001]
         R10 (np.floating):
             Native longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.002]
@@ -32,9 +32,15 @@ def R1_to_C_linear_relaxivity(
             longitudinal relaxation rate, linear with relaxivity model [OSIPI code M.EL1.003]
         - Adapted from equation given in lexicon
     """
-    # Check R1 is a 1D array of floats
-    if not (isinstance(R1, np.ndarray) and R1.ndim == 1 and np.issubdtype(R1.dtype, np.floating)):
-        raise TypeError("R1 must be a 1D NumPy array of np.floating")
+    # Convert input to numpy array with appropriate dtype
+    R1_arr = np.asarray(R1)
+    # Ensure floating-point dtype
+    if not np.issubdtype(R1_arr.dtype, np.floating):
+        R1_arr = R1_arr.astype(np.float64)
+
+    # Check R1 is a 1D array
+    if not (R1_arr.ndim == 1):
+        raise TypeError("R1 must be a 1D array-like object of floating point values")
     elif not (r1 >= 0):
         raise ValueError("r1 must be positive")
-    return (R1 - R10) / r1  # C
+    return (R1_arr - R10) / r1  # C
diff --git a/src/osipi/_signal.py b/src/osipi/_signal.py
index 567737b..3d1e7b7 100644
--- a/src/osipi/_signal.py
+++ b/src/osipi/_signal.py
@@ -1,12 +1,11 @@
 import numpy as np
-from numpy.typing import NDArray
+from numpy.typing import ArrayLike, NDArray
 
 
-# Use a more generic floating-point type annotation
-def signal_linear(R1: NDArray[np.floating], k: np.floating) -> NDArray[np.floating]:
+def signal_linear(R1: ArrayLike, k: np.floating) -> NDArray[np.floating]:
     """Linear model for relationship between R1 and magnitude signal
     Args:
-        R1 (NDArray[np.floating]): longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.001]
+        R1 (ArrayLike): longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.001]
         k (np.floating): proportionality constant in a.u. S [OSIPI code Q.GE1.009]
     Returns:
         NDArray[np.floating]: magnitude signal in a.u. [OSIPI code Q.MS1.001]
@@ -16,20 +15,26 @@ def signal_linear(R1: NDArray[np.floating], k: np.floating) -> NDArray[np.floati
         - OSIPI name: Linear model
         - Adapted from equation given in the Lexicon
     """
+    # Convert input to numpy array with appropriate dtype
+    R1_arr = np.asarray(R1)
+    # Ensure floating-point dtype
+    if not np.issubdtype(R1_arr.dtype, np.floating):
+        R1_arr = R1_arr.astype(np.float64)
+
     # calculate signal
-    return k * R1  # S
+    return k * R1_arr  # S
 
 
 def signal_SPGR(
-    R1: NDArray[np.floating],
-    S0: NDArray[np.floating],
+    R1: ArrayLike,
+    S0: ArrayLike,
     TR: np.floating,
     a: np.floating,
 ) -> NDArray[np.floating]:
     """Steady-state signal for SPGR sequence.
     Args:
-        R1 (NDArray[np.floating]): longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.001]
-        S0 (NDArray[np.floating]): fully T1-relaxed signal in a.u. [OSIPI code Q.MS1.010]
+        R1 (ArrayLike): longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.001]
+        S0 (ArrayLike): fully T1-relaxed signal in a.u. [OSIPI code Q.MS1.010]
         TR (np.floating): repetition time in units of s. [OSIPI code Q.MS1.006]
         a (np.floating): prescribed flip angle in units of deg. [OSIPI code Q.MS1.007]
     Returns:
@@ -40,7 +45,17 @@ def signal_SPGR(
         - OSIPI name: Spoiled gradient recalled echo model
         - Adapted from equation given in the Lexicon and contribution from MJT_UoEdinburgh_UK
     """
+    # Convert inputs to numpy arrays with appropriate dtype
+    R1_arr = np.asarray(R1)
+    S0_arr = np.asarray(S0)
+
+    # Ensure floating-point dtype
+    if not np.issubdtype(R1_arr.dtype, np.floating):
+        R1_arr = R1_arr.astype(np.float64)
+    if not np.issubdtype(S0_arr.dtype, np.floating):
+        S0_arr = S0_arr.astype(np.float64)
+
     # calculate signal
     a_rad = a * np.pi / 180
-    exp_TR_R1 = np.exp(-TR * R1)
-    return S0 * (((1.0 - exp_TR_R1) * np.sin(a_rad)) / (1.0 - exp_TR_R1 * np.cos(a_rad)))  # S
+    exp_TR_R1 = np.exp(-TR * R1_arr)
+    return S0_arr * (((1.0 - exp_TR_R1) * np.sin(a_rad)) / (1.0 - exp_TR_R1 * np.cos(a_rad)))  # S
diff --git a/src/osipi/_signal_to_concentration.py b/src/osipi/_signal_to_concentration.py
index a51dff1..bcaec21 100644
--- a/src/osipi/_signal_to_concentration.py
+++ b/src/osipi/_signal_to_concentration.py
@@ -1,11 +1,11 @@
 import numpy as np
-from numpy.typing import NDArray
+from numpy.typing import ArrayLike, NDArray
 
 from ._electromagnetic_property import R1_to_C_linear_relaxivity
 
 
 def S_to_C_via_R1_SPGR(
-    S: NDArray[np.floating],
+    S: ArrayLike,
     S_baseline: np.floating,
     R10: np.floating,
     TR: np.floating,
@@ -19,7 +19,7 @@ def S_to_C_via_R1_SPGR(
     Converts S -> R1 -> C
 
     Args:
-        S (1D array of np.floating): Vector of magnitude signals in a.u. [OSIPI code Q.MS1.001]
+        S (ArrayLike): Vector of magnitude signals in a.u. [OSIPI code Q.MS1.001]
         S_baseline (np.floating): Pre-contrast magnitude signal in a.u. [OSIPI code Q.MS1.001]
         R10 (np.floating): Native longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.002]
         TR (np.floating): Repetition time in units of s. [OSIPI code Q.MS1.006]
@@ -50,7 +50,7 @@ def S_to_C_via_R1_SPGR(
 
 
 def S_to_R1_SPGR(
-    S: NDArray[np.floating],
+    S: ArrayLike,
     S_baseline: np.floating,
     R10: np.floating,
     TR: np.floating,
@@ -62,7 +62,7 @@ def S_to_R1_SPGR(
     Converts Signal to R1
 
     Args:
-        S (1D array of np.floating): Vector of magnitude signals in a.u. [OSIPI code Q.MS1.001]
+        S (ArrayLike): Vector of magnitude signals in a.u. [OSIPI code Q.MS1.001]
         S_baseline (np.floating): Pre-contrast magnitude signal in a.u. [OSIPI code Q.MS1.001]
         R10 (np.floating): Native longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.002]
         TR (np.floating): Repetition time in units of s. [OSIPI code Q.MS1.006]
@@ -79,9 +79,15 @@ def S_to_R1_SPGR(
           - Forward model: Spoiled gradient recalled echo model [OSIPI code M.SM2.002]
         - Adapted from contribution of LEK_UoEdinburgh_UK
     """
+    # Convert input to numpy array with appropriate dtype
+    S_arr = np.asarray(S)
+    # Ensure floating-point dtype
+    if not np.issubdtype(S_arr.dtype, np.floating):
+        S_arr = S_arr.astype(np.float64)
+
     # Check S is a 1D array of floats
-    if not (isinstance(S, np.ndarray) and S.ndim == 1 and np.issubdtype(S.dtype, np.floating)):
-        raise TypeError("S must be a 1D NumPy array of np.floating")
+    if not (S_arr.ndim == 1):
+        raise TypeError("S must be a 1D array-like object of floating point values")
 
     a_rad = a * np.pi / 180
     # Estimate fully T1-relaxed signal S0 in units of a.u. [OSIPI code Q.MS1.010], then R1
@@ -89,4 +95,4 @@ def S_to_R1_SPGR(
     sin_a = np.sin(a_rad)
     cos_a = np.cos(a_rad)
     S0 = S_baseline * (1 - cos_a * exp_TR_R10) / (sin_a * (1 - exp_TR_R10))
-    return np.log(((S0 * sin_a) - S) / (S0 * sin_a - (S * cos_a))) * (-1 / TR)  # R1
+    return np.log(((S0 * sin_a) - S_arr) / (S0 * sin_a - (S_arr * cos_a))) * (-1 / TR)  # R1
diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py
index 2497539..63dbf43 100755
--- a/src/osipi/_tissue.py
+++ b/src/osipi/_tissue.py
@@ -1,15 +1,15 @@
 import warnings
 
 import numpy as np
-from numpy.typing import NDArray
+from numpy.typing import ArrayLike, NDArray
 from scipy.interpolate import interp1d
 
 from ._convolution import exp_conv
 
 
 def tofts(
-    t: NDArray[np.floating],
-    ca: NDArray[np.floating],
+    t: ArrayLike,
+    ca: ArrayLike,
     Ktrans: np.floating,
     ve: np.floating,
     Ta: np.floating = 30.0,
@@ -18,8 +18,8 @@ def tofts(
     """Tofts model as defined by Tofts and Kermode (1991)
 
     Args:
-        t (NDArray[np.floating]): array of time points in units of sec. [OSIPI code Q.GE1.004]
-        ca (NDArray[np.floating]):
+        t (ArrayLike): array of time points in units of sec. [OSIPI code Q.GE1.004]
+        ca (ArrayLike):
             Arterial concentrations in mM for each time point in t. [OSIPI code Q.IC1.001]
         Ktrans (np.floating):
             Volume transfer constant in units of 1/min. [OSIPI code Q.PH1.008]
@@ -74,14 +74,24 @@ def tofts(
         >>> plt.plot(t, ca, "r", t, ct, "b")
 
     """
-    if not np.allclose(np.diff(t), np.diff(t)[0]):
+    # Convert inputs to numpy arrays with appropriate dtype
+    t_arr = np.asarray(t)
+    ca_arr = np.asarray(ca)
+
+    # Ensure floating-point dtype
+    if not np.issubdtype(t_arr.dtype, np.floating):
+        t_arr = t_arr.astype(np.float64)
+    if not np.issubdtype(ca_arr.dtype, np.floating):
+        ca_arr = ca_arr.astype(np.float64)
+
+    if not np.allclose(np.diff(t_arr), np.diff(t_arr)[0]):
         warnings.warn(
-            ("Non-uniform time spacing detected. Time array may be" " resampled."),
+            ("Non-uniform time spacing detected. Time array may be resampled."),
             stacklevel=2,
         )
 
     if Ktrans <= 0 or ve <= 0:
-        ct = 0 * ca
+        ct = 0 * ca_arr
 
     else:
         # Convert units
@@ -91,54 +101,54 @@ def tofts(
             # Shift the AIF by the arterial delay time (if not zero)
             if Ta != 0:
                 f = interp1d(
-                    t,
-                    ca,
+                    t_arr,
+                    ca_arr,
                     kind="linear",
                     bounds_error=False,
                     fill_value=0,
                 )
-                ca = (t > Ta) * f(t - Ta)
+                ca_arr = (t_arr > Ta) * f(t_arr - Ta)
 
             Tc = ve / Ktrans
-            ct = ve * exp_conv(Tc, t, ca)
+            ct = ve * exp_conv(Tc, t_arr, ca_arr)
 
         else:  # Use convolution by default
             # Calculate the impulse response function
             kep = Ktrans / ve
-            imp = Ktrans * np.exp(-1 * kep * t)
+            imp = Ktrans * np.exp(-1 * kep * t_arr)
 
             # Shift the AIF by the arterial delay time (if not zero)
             if Ta != 0:
                 f = interp1d(
-                    t,
-                    ca,
+                    t_arr,
+                    ca_arr,
                     kind="linear",
                     bounds_error=False,
                     fill_value=0,
                 )
-                ca = (t > Ta) * f(t - Ta)
+                ca_arr = (t_arr > Ta) * f(t_arr - Ta)
 
             # Check if time data grid is uniformly spaced
-            if np.allclose(np.diff(t), np.diff(t)[0]):
+            if np.allclose(np.diff(t_arr), np.diff(t_arr)[0]):
                 # Convolve impulse response with AIF
-                convolution = np.convolve(ca, imp)
+                convolution = np.convolve(ca_arr, imp)
 
                 # Discard unwanted points and make sure time spacing
                 # is correct
-                ct = convolution[0 : len(t)] * t[1]
+                ct = convolution[0 : len(t_arr)] * t_arr[1]
             else:
                 # Resample at the smallest spacing
-                dt = np.min(np.diff(t))
-                t_resampled = np.linspace(t[0], t[-1], int((t[-1] - t[0]) / dt))
+                dt = np.min(np.diff(t_arr))
+                t_resampled = np.linspace(t_arr[0], t_arr[-1], int((t_arr[-1] - t_arr[0]) / dt))
                 ca_func = interp1d(
-                    t,
-                    ca,
+                    t_arr,
+                    ca_arr,
                     kind="quadratic",
                     bounds_error=False,
                     fill_value=0,
                 )
                 imp_func = interp1d(
-                    t,
+                    t_arr,
                     imp,
                     kind="quadratic",
                     bounds_error=False,
@@ -149,8 +159,7 @@ def tofts(
                 # Convolve impulse response with AIF
                 convolution = np.convolve(ca_resampled, imp_resampled)
 
-                # Discard unwanted points and make sure time spacing
-                # is correct
+                # Discard unwanted points and make sure time spacing is correct
                 ct_resampled = convolution[0 : len(t_resampled)] * t_resampled[1]
 
                 # Restore time grid spacing
@@ -161,14 +170,14 @@ def tofts(
                     bounds_error=False,
                     fill_value=0,
                 )
-                ct = ct_func(t)
+                ct = ct_func(t_arr)
 
     return ct
 
 
 def extended_tofts(
-    t: NDArray[np.floating],
-    ca: NDArray[np.floating],
+    t: ArrayLike,
+    ca: ArrayLike,
     Ktrans: np.floating,
     ve: np.floating,
     vp: np.floating,
@@ -178,15 +187,9 @@ def extended_tofts(
     """Extended tofts model as defined by Tofts (1997)
 
     Args:
-        t (NDArray[np.floating]):
+        t (ArrayLike):
             array of time points in units of sec. [OSIPI code Q.GE1.004]
-        ca (NDArray[np.floating]):
-            Arterial concentrations in mM for each time point in t. [OSIPI code Q.IC1.001]
-        Ktrans (np.floating):
-            Volume transfer constant in units of 1/min. [OSIPI code Q.PH1.008]
-        ve (np.floating):
-            Relative volume fraction of the extracellular
-            extravascular compartment (e). [OSIPI code Q.PH1.001.[e]]
+        ca (ArrayLike):
         vp (np.floating):
             Relative volyme fraction of the plasma compartment (p). [OSIPI code Q.PH1.001.[p]]
         Ta (np.floating, optional):
@@ -238,15 +241,24 @@ def extended_tofts(
         >>> plt.plot(t, ca, "r", t, ct, "b")
 
     """
+    # Convert inputs to numpy arrays with appropriate dtype
+    t_arr = np.asarray(t)
+    ca_arr = np.asarray(ca)
+
+    # Ensure floating-point dtype
+    if not np.issubdtype(t_arr.dtype, np.floating):
+        t_arr = t_arr.astype(np.float64)
+    if not np.issubdtype(ca_arr.dtype, np.floating):
+        ca_arr = ca_arr.astype(np.float64)
 
-    if not np.allclose(np.diff(t), np.diff(t)[0]):
+    if not np.allclose(np.diff(t_arr), np.diff(t_arr)[0]):
         warnings.warn(
-            ("Non-uniform time spacing detected. Time array may be" " resampled."),
+            ("Non-uniform time spacing detected. Time array may be resampled."),
             stacklevel=2,
         )
 
     if Ktrans <= 0 or ve <= 0:
-        ct = vp * ca
+        ct = vp * ca_arr
 
     else:
         # Convert units
@@ -256,56 +268,56 @@ def extended_tofts(
             # Shift the AIF by the arterial delay time (if not zero)
             if Ta != 0:
                 f = interp1d(
-                    t,
-                    ca,
+                    t_arr,
+                    ca_arr,
                     kind="linear",
                     bounds_error=False,
                     fill_value=0,
                 )
-                ca = (t > Ta) * f(t - Ta)
+                ca_arr = (t_arr > Ta) * f(t_arr - Ta)
 
             Tc = ve / Ktrans
             # expconv calculates convolution of ca and
             # (1/Tc)exp(-t/Tc), add vp*ca term for extended model
-            ct = (vp * ca) + ve * exp_conv(Tc, t, ca)
+            ct = (vp * ca_arr) + ve * exp_conv(Tc, t_arr, ca_arr)
 
         else:  # Use convolution by default
             # Calculate the impulse response function
             kep = Ktrans / ve
-            imp = Ktrans * np.exp(-1 * kep * t)
+            imp = Ktrans * np.exp(-1 * kep * t_arr)
 
             # Shift the AIF by the arterial delay time (if not zero)
             if Ta != 0:
                 f = interp1d(
-                    t,
-                    ca,
+                    t_arr,
+                    ca_arr,
                     kind="linear",
                     bounds_error=False,
                     fill_value=0,
                 )
-                ca = (t > Ta) * f(t - Ta)
+                ca_arr = (t_arr > Ta) * f(t_arr - Ta)
 
             # Check if time data grid is uniformly spaced
-            if np.allclose(np.diff(t), np.diff(t)[0]):
+            if np.allclose(np.diff(t_arr), np.diff(t_arr)[0]):
                 # Convolve impulse response with AIF
-                convolution = np.convolve(ca, imp)
+                convolution = np.convolve(ca_arr, imp)
 
                 # Discard unwanted points, make sure time spacing is
                 # correct and add vp*ca term for extended model
-                ct = convolution[0 : len(t)] * t[1] + (vp * ca)
+                ct = convolution[0 : len(t_arr)] * t_arr[1] + (vp * ca_arr)
             else:
                 # Resample at the smallest spacing
-                dt = np.min(np.diff(t))
-                t_resampled = np.linspace(t[0], t[-1], int((t[-1] - t[0]) / dt))
+                dt = np.min(np.diff(t_arr))
+                t_resampled = np.linspace(t_arr[0], t_arr[-1], int((t_arr[-1] - t_arr[0]) / dt))
                 ca_func = interp1d(
-                    t,
-                    ca,
+                    t_arr,
+                    ca_arr,
                     kind="quadratic",
                     bounds_error=False,
                     fill_value=0,
                 )
                 imp_func = interp1d(
-                    t,
+                    t_arr,
                     imp,
                     kind="quadratic",
                     bounds_error=False,
@@ -330,6 +342,6 @@ def extended_tofts(
                     bounds_error=False,
                     fill_value=0,
                 )
-                ct = ct_func(t)
+                ct = ct_func(t_arr)
 
     return ct
diff --git a/tests/test_array_like_inputs.py b/tests/test_array_like_inputs.py
new file mode 100644
index 0000000..3a15b59
--- /dev/null
+++ b/tests/test_array_like_inputs.py
@@ -0,0 +1,273 @@
+import numpy as np
+import pytest
+
+import osipi
+
+
+def test_signal_linear_array_like_inputs():
+    # Test with different array-like input types
+
+    # NumPy array input
+    R1 = np.array([0.0, 1.5, 3.0, 4.0, 10.0], dtype=np.float64)
+    k = np.float64(150.0)
+    S_truth = np.array([0.0, 225.0, 450.0, 600.0, 1500.0])
+    S = osipi.signal_linear(R1, k)
+    np.testing.assert_allclose(S_truth, S, rtol=0, atol=1e-7)
+
+    # Python list input
+    R1_list = [0.0, 1.5, 3.0, 4.0, 10.0]
+    S = osipi.signal_linear(R1_list, k)
+    np.testing.assert_allclose(S_truth, S, rtol=0, atol=1e-7)
+
+    # Python tuple input
+    R1_tuple = (0.0, 1.5, 3.0, 4.0, 10.0)
+    S = osipi.signal_linear(R1_tuple, k)
+    np.testing.assert_allclose(S_truth, S, rtol=0, atol=1e-7)
+
+    # Test with float32 dtype
+    R1_32 = np.array([0.0, 1.5, 3.0, 4.0, 10.0], dtype=np.float32)
+    k_32 = np.float32(150.0)
+    S = osipi.signal_linear(R1_32, k_32)
+    np.testing.assert_allclose(S_truth, S, rtol=1e-6, atol=1e-5)  # Reduced precision for float32
+
+    # Test with invalid input
+    with pytest.raises(ValueError):
+        osipi.signal_linear("invalid input", k)
+
+
+def test_signal_SPGR_array_like_inputs():
+    # Test with different array-like input types
+    R1 = np.array([0.1, 0.2, 0.5, 1.0], dtype=np.float64)
+    S0 = np.float64(100)
+    TR = np.float64(5e-3)
+    a = np.float64(15)
+
+    # Calculate expected output with numpy arrays
+    S_truth = osipi.signal_SPGR(R1, S0, TR, a)
+
+    # Python list input for R1
+    R1_list = [0.1, 0.2, 0.5, 1.0]
+    S = osipi.signal_SPGR(R1_list, S0, TR, a)
+    np.testing.assert_allclose(S_truth, S, rtol=0, atol=1e-7)
+
+    # Python tuple input for R1
+    R1_tuple = (0.1, 0.2, 0.5, 1.0)
+    S = osipi.signal_SPGR(R1_tuple, S0, TR, a)
+    np.testing.assert_allclose(S_truth, S, rtol=0, atol=1e-7)
+
+    # Test with array-like S0
+    S0_list = [100, 100, 100, 100]
+    S = osipi.signal_SPGR(R1, S0_list, TR, a)
+    np.testing.assert_allclose(S_truth, S, rtol=0, atol=1e-7)
+
+    # Test with float32 dtype
+    R1_32 = np.array([0.1, 0.2, 0.5, 1.0], dtype=np.float32)
+    S0_32 = np.float32(100)
+    TR_32 = np.float32(5e-3)
+    a_32 = np.float32(15)
+    S = osipi.signal_SPGR(R1_32, S0_32, TR_32, a_32)
+    # Increase tolerance for float32 precision
+    np.testing.assert_allclose(S_truth, S, rtol=1e-4, atol=1e-4)
+
+    # Test with invalid input
+    with pytest.raises(ValueError):
+        osipi.signal_SPGR("invalid input", S0, TR, a)
+
+
+def test_S_to_R1_SPGR_array_like_inputs():
+    # Test data adapted from OSIPI repository tests
+    S_array = np.array([7, 9, 6, 10], dtype=np.float64)
+    S_baseline = np.float64(7)
+    R10 = np.float64(1 / 1.4)
+    TR = np.float64(0.002)
+    a = np.float64(13)
+
+    # Calculate expected result with numpy array
+    R1_truth = osipi.S_to_R1_SPGR(S_array, S_baseline, R10, TR, a)
+
+    # Test with Python list
+    S_list = [7.0, 9.0, 6.0, 10.0]
+    R1 = osipi.S_to_R1_SPGR(S_list, S_baseline, R10, TR, a)
+    np.testing.assert_allclose(R1_truth, R1, rtol=0, atol=1e-7)
+
+    # Test with Python tuple
+    S_tuple = (7.0, 9.0, 6.0, 10.0)
+    R1 = osipi.S_to_R1_SPGR(S_tuple, S_baseline, R10, TR, a)
+    np.testing.assert_allclose(R1_truth, R1, rtol=0, atol=1e-7)
+
+    # Test with float32 dtype
+    S_32 = np.array([7, 9, 6, 10], dtype=np.float32)
+    S_baseline_32 = np.float32(7)
+    R10_32 = np.float32(1 / 1.4)
+    TR_32 = np.float32(0.002)
+    a_32 = np.float32(13)
+    R1 = osipi.S_to_R1_SPGR(S_32, S_baseline_32, R10_32, TR_32, a_32)
+    # Increase tolerance for float32 precision
+    np.testing.assert_allclose(R1_truth, R1, rtol=1e-4, atol=1e-4)
+
+    # Test with non-array input (should raise error as function expects 1D array)
+    with pytest.raises(TypeError):
+        osipi.S_to_R1_SPGR(42, S_baseline, R10, TR, a)
+
+    # Test with 2D array (should raise error as function expects 1D array)
+    with pytest.raises(TypeError):
+        osipi.S_to_R1_SPGR(np.array([[1, 2], [3, 4]]), S_baseline, R10, TR, a)
+
+
+def test_R1_to_C_linear_relaxivity_array_like_inputs():
+    # Test with different array-like input types
+
+    # NumPy array input
+    R1_array = np.array([1, 2, 3, 4, 5, 6], dtype=np.float64)
+    R10 = np.float64(1)
+    r1 = np.float64(5)
+    C_truth = np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], dtype=np.float64)
+    C = osipi.R1_to_C_linear_relaxivity(R1_array, R10, r1)
+    np.testing.assert_allclose(C_truth, C, rtol=0, atol=1e-7)
+
+    # Python list input
+    R1_list = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
+    C = osipi.R1_to_C_linear_relaxivity(R1_list, R10, r1)
+    np.testing.assert_allclose(C_truth, C, rtol=0, atol=1e-7)
+
+    # Python tuple input
+    R1_tuple = (1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
+    C = osipi.R1_to_C_linear_relaxivity(R1_tuple, R10, r1)
+    np.testing.assert_allclose(C_truth, C, rtol=0, atol=1e-7)
+
+    # Test with float32 dtype
+    R1_32 = np.array([1, 2, 3, 4, 5, 6], dtype=np.float32)
+    R10_32 = np.float32(1)
+    r1_32 = np.float32(5)
+    C = osipi.R1_to_C_linear_relaxivity(R1_32, R10_32, r1_32)
+    np.testing.assert_allclose(C_truth, C, rtol=1e-6, atol=1e-6)
+
+    # Test with invalid input types
+    with pytest.raises(TypeError):
+        osipi.R1_to_C_linear_relaxivity(42, R10, r1)  # Not array-like
+
+    with pytest.raises(TypeError):
+        osipi.R1_to_C_linear_relaxivity(np.array([[1, 2], [3, 4]]), R10, r1)  # Not 1D
+
+
+def test_aif_parker_array_like_inputs():
+    # Test with different array-like input types
+
+    # NumPy array input
+    t_array = np.arange(0, 60, 1, dtype=np.float64)
+    ca_array = osipi.aif_parker(t_array)
+
+    # Python list input
+    t_list = list(range(60))
+    ca_list = osipi.aif_parker(t_list)
+    np.testing.assert_allclose(ca_array, ca_list, rtol=1e-7, atol=1e-7)
+
+    # Python tuple input
+    t_tuple = tuple(range(60))
+    ca_tuple = osipi.aif_parker(t_tuple)
+    np.testing.assert_allclose(ca_array, ca_tuple, rtol=1e-7, atol=1e-7)
+
+    # Test with float32 dtype
+    t_32 = np.arange(0, 60, 1, dtype=np.float32)
+    ca_32 = osipi.aif_parker(t_32)
+    np.testing.assert_allclose(ca_array, ca_32, rtol=1e-6, atol=1e-6)
+
+    # Test with invalid input types
+    with pytest.raises(ValueError):
+        osipi.aif_parker("not an array")
+
+
+def test_tofts_array_like_inputs():
+    # Test with different array-like input types
+
+    # NumPy array inputs
+    t_array = np.arange(0, 60, 1, dtype=np.float64)
+    ca_array = osipi.aif_parker(t_array)
+    Ktrans = np.float64(0.6)
+    ve = np.float64(0.2)
+    ct_array = osipi.tofts(t_array, ca_array, Ktrans, ve)
+
+    # Python list inputs
+    t_list = list(range(60))
+    ca_list = list(osipi.aif_parker(t_list))
+    ct_list = osipi.tofts(t_list, ca_list, Ktrans, ve)
+    np.testing.assert_allclose(ct_array, ct_list, rtol=1e-7, atol=1e-7)
+
+    # Python tuple inputs
+    t_tuple = tuple(range(60))
+    ca_tuple = tuple(osipi.aif_parker(t_tuple))
+    ct_tuple = osipi.tofts(t_tuple, ca_tuple, Ktrans, ve)
+    np.testing.assert_allclose(ct_array, ct_tuple, rtol=1e-7, atol=1e-7)
+
+    # Mixed inputs - numpy array and list
+    ct_mixed = osipi.tofts(t_array, ca_list, Ktrans, ve)
+    np.testing.assert_allclose(ct_array, ct_mixed, rtol=1e-7, atol=1e-7)
+
+    # Float32 dtype
+    t_32 = np.arange(0, 60, 1, dtype=np.float32)
+    ca_32 = osipi.aif_parker(t_32)
+    Ktrans_32 = np.float32(0.6)
+    ve_32 = np.float32(0.2)
+    ct_32 = osipi.tofts(t_32, ca_32, Ktrans_32, ve_32)
+    np.testing.assert_allclose(ct_array, ct_32, rtol=1e-5, atol=1e-5)
+
+    # Test with invalid input types
+    with pytest.raises(ValueError):
+        osipi.tofts("invalid", ca_array, Ktrans, ve)
+    with pytest.raises(ValueError):
+        osipi.tofts(t_array, "invalid", Ktrans, ve)
+
+
+def test_extended_tofts_array_like_inputs():
+    # Test with different array-like input types
+
+    # NumPy array inputs
+    t_array = np.arange(0, 60, 1, dtype=np.float64)
+    ca_array = osipi.aif_parker(t_array)
+    Ktrans = np.float64(0.6)
+    ve = np.float64(0.2)
+    vp = np.float64(0.1)
+    ct_array = osipi.extended_tofts(t_array, ca_array, Ktrans, ve, vp)
+
+    # Python list inputs
+    t_list = list(range(60))
+    ca_list = list(osipi.aif_parker(t_list))
+    ct_list = osipi.extended_tofts(t_list, ca_list, Ktrans, ve, vp)
+    np.testing.assert_allclose(ct_array, ct_list, rtol=1e-7, atol=1e-7)
+
+    # Python tuple inputs
+    t_tuple = tuple(range(60))
+    ca_tuple = tuple(osipi.aif_parker(t_tuple))
+    ct_tuple = osipi.extended_tofts(t_tuple, ca_tuple, Ktrans, ve, vp)
+    np.testing.assert_allclose(ct_array, ct_tuple, rtol=1e-7, atol=1e-7)
+
+    # Mixed inputs - numpy array and list
+    ct_mixed = osipi.extended_tofts(t_array, ca_list, Ktrans, ve, vp)
+    np.testing.assert_allclose(ct_array, ct_mixed, rtol=1e-7, atol=1e-7)
+
+    # Float32 dtype
+    t_32 = np.arange(0, 60, 1, dtype=np.float32)
+    ca_32 = osipi.aif_parker(t_32)
+    Ktrans_32 = np.float32(0.6)
+    ve_32 = np.float32(0.2)
+    vp_32 = np.float32(0.1)
+    ct_32 = osipi.extended_tofts(t_32, ca_32, Ktrans_32, ve_32, vp_32)
+    np.testing.assert_allclose(ct_array, ct_32, rtol=1e-5, atol=1e-5)
+
+    # Test with invalid input types
+    with pytest.raises(ValueError):
+        osipi.extended_tofts("invalid", ca_array, Ktrans, ve, vp)
+    with pytest.raises(ValueError):
+        osipi.extended_tofts(t_array, "invalid", Ktrans, ve, vp)
+
+
+if __name__ == "__main__":
+    test_signal_linear_array_like_inputs()
+    test_signal_SPGR_array_like_inputs()
+    test_S_to_R1_SPGR_array_like_inputs()
+    test_R1_to_C_linear_relaxivity_array_like_inputs()
+    test_aif_parker_array_like_inputs()
+    test_tofts_array_like_inputs()
+    test_extended_tofts_array_like_inputs()
+
+    print("All array-like input tests passed!!")