diff --git a/src/sp2xr/apply_calib.py b/src/sp2xr/apply_calib.py
deleted file mode 100644
index 3c23ea4..0000000
--- a/src/sp2xr/apply_calib.py
+++ /dev/null
@@ -1,42 +0,0 @@
-import numpy as np
-import pandas as pd
-from sp2xr.calibration import calibrate_particle_data
-from sp2xr.flag_single_particle_data import define_flags, add_thin_thick_flags
-
-
-def calibrate_single_particle(ddf_raw, instr_config, run_config):
-    meta_cal = ddf_raw._meta.copy()
-    meta_cal["BC mass"] = pd.Series([], dtype="float64")
-    meta_cal["Opt diam"] = pd.Series([], dtype="float64")
-    meta_cal["time_lag"] = pd.Series([], dtype="int8")
-    meta_cal["date"] = pd.Series([], dtype="datetime64[ns]")
-    meta_cal["hour"] = pd.Series([], dtype="int8")
-
-    ddf_cal = ddf_raw.map_partitions(
-        calibrate_particle_data, config=run_config, meta=meta_cal
-    )
-    ddf_cal["time_lag"] = ddf_cal["Incand Peak Time"] - ddf_cal["Scatter Peak Time"]
-    ddf_cal["ratio_inc_scatt"] = np.log10(ddf_cal["Incand relPeak"]) / np.log10(
-        ddf_cal["Scatter relPeak"]
-    )
-
-    ddf_cal = define_flags(ddf_cal, instr_config, run_config)
-
-    meta_cnts = add_thin_thick_flags(ddf_cal._meta)  # builds empty columns for meta
-    ddf_cal = ddf_cal.map_partitions(add_thin_thick_flags, meta=meta_cnts)
-
-    # Create two new columns where the particles out of range are changed to np.nan (!! do not change them to 0, otherwise below where we resample and count, the 0s would be counted being a not null value!!)
-    ddf_cal["BC mass within range"] = ddf_cal["BC mass"].where(
-        ddf_cal["flag_valid_inc_signal_in_range"], np.nan
-    )
-    ddf_cal["Opt diam scatt only within range"] = ddf_cal["Opt diam scatt only"].where(
-        (
-            ddf_cal["flag_valid_scatt_signal_in_range"]
-            & ~ddf_cal["flag_valid_inc_signal"]
-        ),
-        np.nan,
-    )
-
-    ddf_cal["temporary_col"] = 1
-
-    return ddf_cal
diff --git a/src/sp2xr/calibration.py b/src/sp2xr/calibration.py
index d94d147..31b18de 100644
--- a/src/sp2xr/calibration.py
+++ b/src/sp2xr/calibration.py
@@ -1,16 +1,35 @@
+"""
+SP2-XR calibration functions for single particle data processing.
+
+This module provides functions for applying scattering and incandescence calibrations
+to SP2-XR particle data, as well as conversion functions between mass and diameter.
+"""
+
+from __future__ import annotations
+
+from typing import Any
 import numpy as np
-from numpy.polynomial import Polynomial
+import pandas as pd
+import dask.dataframe as dd
+from sp2xr.flag_single_particle_data import define_flags, add_thin_thick_flags
+from .calibration_constants import (
+    SP2_calibCurveSpline,
+    Aquadag_RhoVsLogMspline_Var25,
+    Aquadag_RhoVsLogMspline_Var8,
+    Fullerene_RhoVsLogMspline_Var2,
+    Fullerene_RhoVsLogMspline_Var5,
+    Fullerene_RhoVsLogMspline_Var8,
+    Fullerene_RhoVsLogMspline_Var8_old,
+    Fullerene_RhoVsLogDspline_Var2,
+    Fullerene_RhoVsLogDspline_Var5,
+    Fullerene_RhoVsLogDspline_Var8,
+    Fullerene_RhoVsLogDspline_Var8_old,
+    GlassyCarbonAlpha_Mass2Diam,
+)
 
 
-def apply_calibration(partition, config):
-    """Apply calibration to a pandas partition."""
-    pdf = calibrate_particle_data(partition, config)
-    # pdf["date"] = partition.index.date.astype("datetime64[ns]")
-    # pdf["hour"] = partition.index.hour.astype("int8")
-    return pdf
-
-
-def polynomial(x, *coefficients):
+# Core mathematical functions
+def polynomial(x: float | np.ndarray, *coefficients: float) -> float | np.ndarray:
     """
     Function to generate a generic polynomial curve. The number of input coefficicnets define the degree of the polynomial curve
 
@@ -33,7 +52,7 @@ def polynomial(x, *coefficients):
     return result
 
 
-def powerlaw(x, *coeff):
+def powerlaw(x: float | np.ndarray, *coeff: float) -> float | np.ndarray:
     """
     Generic power-law calibration function.
 
@@ -52,17 +71,22 @@ def powerlaw(x, *coeff):
     return A * (x**B) + C
 
 
-def apply_inc_calibration(df, curve_type=None, params=None):
+# Low-level calibration application functions
+def apply_inc_calibration(
+    df: pd.DataFrame | dd.DataFrame,
+    curve_type: str | None = None,
+    params: list[float] | None = None,
+) -> pd.DataFrame | dd.DataFrame:
     """
-    Apply incandescence calibration to a Dask DataFrame.
+    Apply incandescence calibration to a DataFrame.
 
     Parameters:
-        df (pd.DataFrame): particle dataframe
+        df (pd.DataFrame | dd.DataFrame): particle dataframe
         curve_type (str): 'polynomial', 'powerlaw', or None (use raw mass)
         params (list or tuple): calibration parameters
 
     Returns:
-        df (pd.DataFrame) with new column 'BC mass'
+        df (pd.DataFrame | dd.DataFrame) with new column 'BC mass'
     """
     inc = df["Incand relPeak"].to_numpy()
 
@@ -85,17 +109,21 @@ def apply_inc_calibration(df, curve_type=None, params=None):
     return df
 
 
-def apply_scatt_calibration(df, curve_type=None, params=None):
+def apply_scatt_calibration(
+    df: pd.DataFrame | dd.DataFrame,
+    curve_type: str | None = None,
+    params: list[float] | None = None,
+) -> pd.DataFrame | dd.DataFrame:
     """
-    Apply scattering calibration to a Dask DataFrame.
+    Apply scattering calibration to a DataFrame.
 
     Parameters:
-        df (pd.DataFrame)
+        df (pd.DataFrame | dd.DataFrame): particle dataframe
         curve_type (str): 'polynomial', 'powerlaw', or None (use raw size)
         params (list or tuple): calibration parameters
 
     Returns:
-        df (pd.DataFrame) with new column 'Opt diam'
+        df (pd.DataFrame | dd.DataFrame) with new column 'Opt diam'
     """
     scatt = df["Scatter relPeak"].to_numpy()
 
@@ -118,8 +146,27 @@ def apply_scatt_calibration(df, curve_type=None, params=None):
     return df
 
 
-def calibrate_particle_data(df, config):
-    """Wrapper to apply both incandescence and scattering calibration based on config."""
+# Mid-level calibration wrapper
+def calibrate_particle_data(
+    df: pd.DataFrame | dd.DataFrame, config: dict[str, Any]
+) -> pd.DataFrame | dd.DataFrame:
+    """
+    Apply both incandescence and scattering calibration based on config.
+
+    This function works with both pandas and dask DataFrames.
+
+    Parameters
+    ----------
+    df : pd.DataFrame | dd.DataFrame
+        Particle dataframe
+    config : dict[str, Any]
+        Configuration dictionary containing calibration parameters
+
+    Returns
+    -------
+    pd.DataFrame | dd.DataFrame
+        Calibrated dataframe with 'BC mass', 'Opt diam', 'date', and 'hour' columns
+    """
     cal_cfg = config.get("calibration", {})
     inc_conf = cal_cfg.get("incandescence", {})
     scatt_conf = cal_cfg.get("scattering", {})
@@ -140,10 +187,82 @@ def calibrate_particle_data(df, config):
     return df
 
 
-# %% Convertions mass to D , N to mass
+# High-level calibration function (merged from apply_calib.py)
+def calibrate_single_particle(
+    ddf_raw: dd.DataFrame, instr_config: dict[str, Any], run_config: dict[str, Any]
+) -> dd.DataFrame:
+    """
+    Apply calibration to single particle data with flags and metadata processing.
+
+    This function performs the complete calibration workflow:
+    1. Apply scattering and incandescence calibrations
+    2. Calculate derived quantities (time_lag, ratio_inc_scatt)
+    3. Apply quality control flags
+    4. Add thin/thick flags
+    5. Create range-filtered columns
+
+    Parameters
+    ----------
+    ddf_raw : dask.DataFrame
+        Raw particle data
+    instr_config : dict
+        Instrument configuration
+    run_config : dict
+        Run configuration including calibration parameters
+
+    Returns
+    -------
+    ddf_cal : dask.DataFrame
+        Calibrated particle data with flags and derived quantities
+    """
+    meta_cal = ddf_raw._meta.copy()
+    meta_cal["BC mass"] = pd.Series([], dtype="float64")
+    meta_cal["Opt diam"] = pd.Series([], dtype="float64")
+    meta_cal["time_lag"] = pd.Series([], dtype="int8")
+    meta_cal["date"] = pd.Series([], dtype="datetime64[ns]")
+    meta_cal["hour"] = pd.Series([], dtype="int8")
+
+    ddf_cal = ddf_raw.map_partitions(
+        calibrate_particle_data, config=run_config, meta=meta_cal
+    )
+    ddf_cal["time_lag"] = ddf_cal["Incand Peak Time"] - ddf_cal["Scatter Peak Time"]
+    ddf_cal["ratio_inc_scatt"] = np.log10(ddf_cal["Incand relPeak"]) / np.log10(
+        ddf_cal["Scatter relPeak"]
+    )
+
+    ddf_cal = define_flags(ddf_cal, instr_config, run_config)
+
+    meta_cnts = add_thin_thick_flags(ddf_cal._meta)  # builds empty columns for meta
+    ddf_cal = ddf_cal.map_partitions(add_thin_thick_flags, meta=meta_cnts)
+
+    # Create two new columns where the particles out of range are changed to np.nan (!! do not change them to 0, otherwise below where we resample and count, the 0s would be counted being a not null value!!)
+    ddf_cal["BC mass within range"] = ddf_cal["BC mass"].where(
+        ddf_cal["flag_valid_inc_signal_in_range"], np.nan
+    )
+    ddf_cal["Opt diam scatt only within range"] = ddf_cal["Opt diam scatt only"].where(
+        (
+            ddf_cal["flag_valid_scatt_signal_in_range"]
+            & ~ddf_cal["flag_valid_inc_signal"]
+        ),
+        np.nan,
+    )
+
+    ddf_cal["temporary_col"] = 1
+
+    return ddf_cal
 
 
-def mass2meqDiam(mass_array, rho=1800):
+# Legacy wrapper for backward compatibility
+def apply_calibration(partition: pd.DataFrame, config: dict[str, Any]) -> pd.DataFrame:
+    """Apply calibration to a pandas partition."""
+    pdf = calibrate_particle_data(partition, config)
+    # pdf["date"] = partition.index.date.astype("datetime64[ns]")
+    # pdf["hour"] = partition.index.hour.astype("int8")
+    return pdf
+
+
+# Conversion functions
+def mass2meqDiam(mass_array: list[float] | np.ndarray, rho: float = 1800) -> np.ndarray:
     """
     Calculate mass equivalent diameters for the input array of mass values.
     inputs: mass_array values in fg, particle density in kg m-3 (rho)
@@ -153,7 +272,9 @@ def mass2meqDiam(mass_array, rho=1800):
     return Dp
 
 
-def dNdlogDp_to_dMdlogDp(N, Dp, rho=1800):
+def dNdlogDp_to_dMdlogDp(
+    N: np.ndarray, Dp: np.ndarray, rho: float = 1800
+) -> np.ndarray:
     """This is from functions.py from Rob"""
     """
     Given a normalized number size distribution [dN/dlogDp, cm^-3] defined over the diameters given by Dp [nm]
@@ -166,7 +287,9 @@ def dNdlogDp_to_dMdlogDp(N, Dp, rho=1800):
     return M
 
 
-def BC_mass_to_diam(mass, rho_eff=1800, BC_type=""):
+def BC_mass_to_diam(
+    mass: float | np.ndarray, rho_eff: float = 1800, BC_type: str = ""
+) -> tuple[float | np.ndarray, float | np.ndarray]:
     """
     Thi function calculates the effective density and mobility (volume equivalent?)  diameter
     for BC calibration materials applying different parameterizations.
@@ -231,8 +354,28 @@ def BC_mass_to_diam(mass, rho_eff=1800, BC_type=""):
     return density, diam
 
 
-def BC_diam_to_mass(diam, rho_eff=1800, BC_type=""):
+def BC_diam_to_mass(
+    diam: float | np.ndarray, rho_eff: float = 1800, BC_type: str = ""
+) -> tuple[float | np.ndarray, float | np.ndarray]:
+    """
+    Convert BC diameter to mass using various calibration materials.
 
+    Parameters
+    ----------
+    diam : float or array
+        Diameter in nm
+    rho_eff : float, optional
+        Effective density [kg/m3]. Default is 1800.
+    BC_type : str, optional
+        Calibration material type. Default is "".
+
+    Returns
+    -------
+    density : float or array
+        Effective density [kg/m3]
+    mass : float or array
+        Mass in fg
+    """
     if BC_type == "":  # pers comm Kondo 2009, spline fit
         SplineCoefD2M, nSplineSegmD2M = Fullerene_RhoVsLogDspline_Var2(SplineCoef=[])
         density = SP2_calibCurveSpline(diam, SplineCoefD2M)
@@ -254,335 +397,3 @@ def BC_diam_to_mass(diam, rho_eff=1800, BC_type=""):
     mass = 1e-9 * density * np.pi / 6 * np.array(diam) ** 3
 
     return density, mass
-
-
-def SP2_calibCurveSpline(mass, SplineCoef):
-    """
-    Given the mass and the spline coefficients list, it calculates the corresponding density.
-
-    Parameters
-    ----------
-    mass : float or array of floats
-        BC mass [fg].
-    SplineCoefM2D : list of floats
-        Coefficients and lowercuts of th espline to be used.
-
-    Returns
-    -------
-    values : float or array of floats
-        BC effective density [kg/m3] corresponding to the mass and specific parameterization.
-
-    """
-    coefficients_array = np.array(SplineCoef)
-    segments = np.delete(coefficients_array, np.arange(3, len(coefficients_array), 4))
-    segments_2 = np.array([segments[i : i + 3] for i in range(0, len(segments), 3)])[
-        ::-1
-    ]
-    lowercuts = coefficients_array[3::4]
-    sorted_lowercuts = np.sort(
-        np.append(lowercuts, [np.log10(mass[0]), np.log10(mass[-1])])
-    )
-
-    values = np.zeros_like(mass)  # Initialize an array to store interpolated values
-
-    for i in range(len(sorted_lowercuts) - 1):
-        if i == (len(sorted_lowercuts) - 2):
-            mask = (np.log10(mass) >= sorted_lowercuts[i]) & (
-                np.log10(mass) <= sorted_lowercuts[i + 1]
-            )
-        else:
-            mask = (np.log10(mass) >= sorted_lowercuts[i]) & (
-                np.log10(mass) < sorted_lowercuts[i + 1]
-            )
-        if np.any(mask):
-            poly = Polynomial(segments_2[i])
-            # values[mask] = poly(np.log10(mass[mask]))
-            values[mask] = poly(np.log10(np.array(mass)[mask]))
-
-    return values
-
-
-def Aquadag_RhoVsLogMspline_Var25(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    SplineCoef.append(1392.54)
-    SplineCoef.append(-618.898)
-    SplineCoef.append(104.364)
-    SplineCoef.append(1.87)
-    SplineCoef.append(689)
-    SplineCoef.append(133.555)
-    SplineCoef.append(-96.8268)
-    SplineCoef.append(0.64)
-    SplineCoef.append(785.701)
-    SplineCoef.append(-168.635)
-    SplineCoef.append(139.258)
-    return SplineCoef, int(
-        SP2_calibCurveSplineCheck(SplineCoef)
-    )  # SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Aquadag_RhoVsLogMspline_Var8(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    SplineCoef.append(1662.41)
-    SplineCoef.append(-1087.41)
-    SplineCoef.append(214.636)
-    SplineCoef.append(1.7)
-    SplineCoef.append(689.745)
-    SplineCoef.append(56.897)
-    SplineCoef.append(-121.925)
-    SplineCoef.append(0.33)
-    SplineCoef.append(720.248)
-    SplineCoef.append(-127.967)
-    SplineCoef.append(158.172)
-    SplineCoef.append(-1.23)
-    SplineCoef.append(299.662)
-    SplineCoef.append(-811.846)
-    SplineCoef.append(-119.828)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Fullerene_RhoVsLogMspline_Var2(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    SplineCoef.append(1180.1)
-    SplineCoef.append(-540.043)
-    SplineCoef.append(106.852)
-    SplineCoef.append(1.8)
-    SplineCoef.append(713.102)
-    SplineCoef.append(-21.1519)
-    SplineCoef.append(-37.2841)
-    SplineCoef.append(0.55)
-    SplineCoef.append(916.24)
-    SplineCoef.append(-759.835)
-    SplineCoef.append(634.246)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Fullerene_RhoVsLogMspline_Var5(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    SplineCoef.append(1539.19)
-    SplineCoef.append(-904.484)
-    SplineCoef.append(189.201)
-    SplineCoef.append(1.83)
-    SplineCoef.append(503.999)
-    SplineCoef.append(226.877)
-    SplineCoef.append(-119.914)
-    SplineCoef.append(0.8)
-    SplineCoef.append(718.903)
-    SplineCoef.append(-310.384)
-    SplineCoef.append(215.874)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Fullerene_RhoVsLogMspline_Var8(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    SplineCoef.append(62.7393)
-    SplineCoef.append(781.804)
-    SplineCoef.append(-325.113)
-    SplineCoef.append(1.6)
-    SplineCoef.append(652.935)
-    SplineCoef.append(44.0588)
-    SplineCoef.append(-94.5682)
-    SplineCoef.append(0.6)
-    SplineCoef.append(750.496)
-    SplineCoef.append(-281.145)
-    SplineCoef.append(176.435)
-    SplineCoef.append(-0.7)
-    SplineCoef.append(573.019)
-    SplineCoef.append(-788.222)
-    SplineCoef.append(-185.763)
-    SplineCoef.append(-1.7)
-    SplineCoef.append(1016.31)
-    SplineCoef.append(-266.707)
-    SplineCoef.append(-32.3765)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Fullerene_RhoVsLogMspline_Var8_old(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    SplineCoef.append(984.914)
-    SplineCoef.append(-468.01)
-    SplineCoef.append(101)
-    SplineCoef.append(1.19)
-    SplineCoef.append(568.667)
-    SplineCoef.append(231.565)
-    SplineCoef.append(-192.939)
-    SplineCoef.append(0.715)
-    SplineCoef.append(747.715)
-    SplineCoef.append(-269.27)
-    SplineCoef.append(157.295)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def GlassyCarbonAlpha_Mass2Diam(mass):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    diam = 1e3 * (6 * mass / np.pi / 1420) ** (1 / 3)
-    return diam
-
-
-def Fullerene_RhoVsLogDspline_Var2(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
-    # units: rho [kg/m³]; diameter [nm]
-    #  experimental mobility/mass data taken from: personal communication with the Kondo research group on 04/09/2009;
-    SplineCoef.append(8799.39)
-    SplineCoef.append(-5478.89)
-    SplineCoef.append(904.097)
-    SplineCoef.append(2.8)
-    SplineCoef.append(-115.387)
-    SplineCoef.append(888.81)
-    SplineCoef.append(-232.992)
-    SplineCoef.append(2.35)
-    SplineCoef.append(17587.6)
-    SplineCoef.append(-14177.6)
-    SplineCoef.append(2972.62)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Fullerene_RhoVsLogDspline_Var5(SplineCoef=[]):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
-    # units: rho [kg/m³]; diameter [nm]
-    # note: experimental mobility/mass data taken from: Moteki et al., Aerosol Sci. Technol., 44, 663-675, 2010;
-    SplineCoef.append(13295.3)
-    SplineCoef.append(-8550.29)
-    SplineCoef.append(1423.82)
-    SplineCoef.append(2.8)
-    SplineCoef.append(-4979.54)
-    SplineCoef.append(4503.18)
-    SplineCoef.append(-907.154)
-    SplineCoef.append(2.45)
-    SplineCoef.append(8500)
-    SplineCoef.append(-6500.53)
-    SplineCoef.append(1338.5)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Fullerene_RhoVsLogDspline_Var8(SplineCoef):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
-    # units: rho [kg/m³]; diameter [nm]
-    # note: experimental mobility/mass data taken from: PSI data for fullerene soot sample obtained from Shuka Schwarz, measured at ETH in Oct. 2010;
-    SplineCoef.append(-2281.86)
-    SplineCoef.append(2687.72)
-    SplineCoef.append(-613.973)
-    SplineCoef.append(2.735)
-    SplineCoef.append(-1987.85)
-    SplineCoef.append(2472.72)
-    SplineCoef.append(-574.667)
-    SplineCoef.append(2.35)
-    SplineCoef.append(7319.9)
-    SplineCoef.append(-5448.76)
-    SplineCoef.append(1110.76)
-    SplineCoef.append(1.85)
-    SplineCoef.append(-520.822)
-    SplineCoef.append(3027.69)
-    SplineCoef.append(-1180.18)
-    SplineCoef.append(1.48)
-    SplineCoef.append(1421.03)
-    SplineCoef.append(403.564)
-    SplineCoef.append(-293.649)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Fullerene_RhoVsLogDspline_Var8_old(SplineCoef):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
-    # units: rho [kg/m³]; diameter [nm]
-    # note: experimental mobility/mass data taken from: PSI data for fullerene soot sample obtained from Shuka Schwarz, measured at ETH in Oct. 2010;
-    SplineCoef.append(6744.81)
-    SplineCoef.append(-4200.82)
-    SplineCoef.append(700)
-    SplineCoef.append(2.57)
-    SplineCoef.append(-6902.08)
-    SplineCoef.append(6419.32)
-    SplineCoef.append(-1366.18)
-    SplineCoef.append(2.4)
-    SplineCoef.append(7288.07)
-    SplineCoef.append(-5405.8)
-    SplineCoef.append(1097.39)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Aquadag_RhoVsLogDspline_Var25(SplineCoef):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
-    # units: rho [kg/m³]; diameter [nm]
-    # note: experimental mobility/mass data from: Moteki et al., Aerosol Sci. Technol., 44, 663-675, 2010 (and they agree with pers. comm. from 04/09/2009)
-    SplineCoef.append(10000)
-    SplineCoef.append(-6088.15)
-    SplineCoef.append(974.717)
-    SplineCoef.append(2.8)
-    SplineCoef.append(-3200)
-    SplineCoef.append(3340.42)
-    SplineCoef.append(-708.956)
-    SplineCoef.append(2.35)
-    SplineCoef.append(6490.76)
-    SplineCoef.append(-4907.03)
-    SplineCoef.append(1045.82)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Aquadag_RhoVsLogDspline_Var8(SplineCoef):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
-    # units: rho [kg/m³]; diameter [nm]
-    # note: experimental mobility/mass data taken from: "grand average" of PSI's (APM at ETH in Oct. 2010) and DMT's (Subu; CPMA in 2011 by XXXX) data for classic and new aquadag batches;
-    SplineCoef.append(10092.9)
-    SplineCoef.append(-6171.83)
-    SplineCoef.append(970.172)
-    SplineCoef.append(2.75)
-    SplineCoef.append(-2627.45)
-    SplineCoef.append(3079.33)
-    SplineCoef.append(-711.856)
-    SplineCoef.append(2.25)
-    SplineCoef.append(6140.89)
-    SplineCoef.append(-4714.75)
-    SplineCoef.append(1020.16)
-    SplineCoef.append(1.65)
-    SplineCoef.append(1063.19)
-    SplineCoef.append(1440.04)
-    SplineCoef.append(-844.928)
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def Aquadag_RhoVsLogDspline_Var8_old(SplineCoef):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
-    # units: rho [kg/m³]; diameter [nm]
-    # note: experimental mobility/mass data taken from: "grand average" of PSI's (APM at ETH in Oct. 2010) and DMT's (Subu; CPMA in 2011 by XXXX) data for classic and new aquadag batches;
-    SplineCoef[0] = 7900.57
-    SplineCoef[1] = -4964.2
-    SplineCoef[2] = 830.317
-    SplineCoef[3] = 2.6
-    SplineCoef[4] = -5371.15
-    SplineCoef[5] = 5244.82
-    SplineCoef[6] = -1132.96
-    SplineCoef[7] = 2.35
-    SplineCoef[8] = 5809.57
-    SplineCoef[9] = -4270.69
-    SplineCoef[10] = 891.62
-    SplineCoef[11] = 1.6
-    SplineCoef[12] = 904.355
-    SplineCoef[13] = 1860.83
-    SplineCoef[14] = -1024.48
-    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
-
-
-def GlassyCarbonAlpha_Diam2Mass(DiamMob):
-    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
-    # This function calculates the mass of a glassy carbon particle of a given mobility diameter.
-    # exact description: Glassy carbon spherical powder, 0.4-12 micron, type 2, Alpha Aesar Art. No. 38008
-    # note: the bulk density is taken from the supplier's data sheet and an email by Ronald Becht (Technical service, Alfa Aesar)
-    kGlassyCarbonAlphaRho = 1420  # kg/m3
-    mass = kGlassyCarbonAlphaRho * np.pi / 6 * DiamMob**3 * 1e-9
-    return mass
-
-
-def SP2_calibCurveSplineCheck(CalibCoef):
-    # this function comes from the sp2 toolkit
-    # This function checks the dimension of the calibration coefficient wave.
-    nCoef = len(CalibCoef)
-    if np.mod(nCoef, 4) != 3:
-        raise ValueError(
-            "The number of coefficients handed over to the function must be 3,7,11, or ...!"
-        )
-    return np.ceil(nCoef / 4)
diff --git a/src/sp2xr/calibration_constants.py b/src/sp2xr/calibration_constants.py
new file mode 100644
index 0000000..9559eac
--- /dev/null
+++ b/src/sp2xr/calibration_constants.py
@@ -0,0 +1,400 @@
+"""
+Calibration constants and spline coefficient functions for SP2-XR data processing.
+
+This module contains spline coefficient functions for different BC calibration materials
+and their parameterizations. These functions come from the original SP2 toolkit.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+from numpy.polynomial import Polynomial
+
+
+def SP2_calibCurveSpline(
+    mass: float | np.ndarray, SplineCoef: list[float]
+) -> np.ndarray:
+    """
+    Given the mass and the spline coefficients list, it calculates the corresponding density.
+
+    Parameters
+    ----------
+    mass : float or array of floats
+        BC mass [fg].
+    SplineCoefM2D : list of floats
+        Coefficients and lowercuts of th espline to be used.
+
+    Returns
+    -------
+    values : float or array of floats
+        BC effective density [kg/m3] corresponding to the mass and specific parameterization.
+
+    """
+    coefficients_array = np.array(SplineCoef)
+    segments = np.delete(coefficients_array, np.arange(3, len(coefficients_array), 4))
+    segments_2 = np.array([segments[i : i + 3] for i in range(0, len(segments), 3)])[
+        ::-1
+    ]
+    lowercuts = coefficients_array[3::4]
+    sorted_lowercuts = np.sort(
+        np.append(lowercuts, [np.log10(mass[0]), np.log10(mass[-1])])
+    )
+
+    values = np.zeros_like(mass)  # Initialize an array to store interpolated values
+
+    for i in range(len(sorted_lowercuts) - 1):
+        if i == (len(sorted_lowercuts) - 2):
+            mask = (np.log10(mass) >= sorted_lowercuts[i]) & (
+                np.log10(mass) <= sorted_lowercuts[i + 1]
+            )
+        else:
+            mask = (np.log10(mass) >= sorted_lowercuts[i]) & (
+                np.log10(mass) < sorted_lowercuts[i + 1]
+            )
+        if np.any(mask):
+            poly = Polynomial(segments_2[i])
+            # values[mask] = poly(np.log10(mass[mask]))
+            values[mask] = poly(np.log10(np.array(mass)[mask]))
+
+    return values
+
+
+def SP2_calibCurveSplineCheck(CalibCoef: list[float]) -> int:
+    # this function comes from the sp2 toolkit
+    # This function checks the dimension of the calibration coefficient wave.
+    nCoef = len(CalibCoef)
+    if np.mod(nCoef, 4) != 3:
+        raise ValueError(
+            "The number of coefficients handed over to the function must be 3,7,11, or ...!"
+        )
+    return int(np.ceil(nCoef / 4))
+
+
+# Mass-to-density spline functions
+def Aquadag_RhoVsLogMspline_Var25(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(1392.54)
+    SplineCoef.append(-618.898)
+    SplineCoef.append(104.364)
+    SplineCoef.append(1.87)
+    SplineCoef.append(689)
+    SplineCoef.append(133.555)
+    SplineCoef.append(-96.8268)
+    SplineCoef.append(0.64)
+    SplineCoef.append(785.701)
+    SplineCoef.append(-168.635)
+    SplineCoef.append(139.258)
+    return SplineCoef, int(
+        SP2_calibCurveSplineCheck(SplineCoef)
+    )  # SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Aquadag_RhoVsLogMspline_Var8(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(1662.41)
+    SplineCoef.append(-1087.41)
+    SplineCoef.append(214.636)
+    SplineCoef.append(1.7)
+    SplineCoef.append(689.745)
+    SplineCoef.append(56.897)
+    SplineCoef.append(-121.925)
+    SplineCoef.append(0.33)
+    SplineCoef.append(720.248)
+    SplineCoef.append(-127.967)
+    SplineCoef.append(158.172)
+    SplineCoef.append(-1.23)
+    SplineCoef.append(299.662)
+    SplineCoef.append(-811.846)
+    SplineCoef.append(-119.828)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Fullerene_RhoVsLogMspline_Var2(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(1180.1)
+    SplineCoef.append(-540.043)
+    SplineCoef.append(106.852)
+    SplineCoef.append(1.8)
+    SplineCoef.append(713.102)
+    SplineCoef.append(-21.1519)
+    SplineCoef.append(-37.2841)
+    SplineCoef.append(0.55)
+    SplineCoef.append(916.24)
+    SplineCoef.append(-759.835)
+    SplineCoef.append(634.246)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Fullerene_RhoVsLogMspline_Var5(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(1539.19)
+    SplineCoef.append(-904.484)
+    SplineCoef.append(189.201)
+    SplineCoef.append(1.83)
+    SplineCoef.append(503.999)
+    SplineCoef.append(226.877)
+    SplineCoef.append(-119.914)
+    SplineCoef.append(0.8)
+    SplineCoef.append(718.903)
+    SplineCoef.append(-310.384)
+    SplineCoef.append(215.874)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Fullerene_RhoVsLogMspline_Var8(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(62.7393)
+    SplineCoef.append(781.804)
+    SplineCoef.append(-325.113)
+    SplineCoef.append(1.6)
+    SplineCoef.append(652.935)
+    SplineCoef.append(44.0588)
+    SplineCoef.append(-94.5682)
+    SplineCoef.append(0.6)
+    SplineCoef.append(750.496)
+    SplineCoef.append(-281.145)
+    SplineCoef.append(176.435)
+    SplineCoef.append(-0.7)
+    SplineCoef.append(573.019)
+    SplineCoef.append(-788.222)
+    SplineCoef.append(-185.763)
+    SplineCoef.append(-1.7)
+    SplineCoef.append(1016.31)
+    SplineCoef.append(-266.707)
+    SplineCoef.append(-32.3765)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Fullerene_RhoVsLogMspline_Var8_old(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(984.914)
+    SplineCoef.append(-468.01)
+    SplineCoef.append(101)
+    SplineCoef.append(1.19)
+    SplineCoef.append(568.667)
+    SplineCoef.append(231.565)
+    SplineCoef.append(-192.939)
+    SplineCoef.append(0.715)
+    SplineCoef.append(747.715)
+    SplineCoef.append(-269.27)
+    SplineCoef.append(157.295)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+# Diameter-to-density spline functions
+def Fullerene_RhoVsLogDspline_Var2(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
+    # units: rho [kg/m³]; diameter [nm]
+    #  experimental mobility/mass data taken from: personal communication with the Kondo research group on 04/09/2009;
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(8799.39)
+    SplineCoef.append(-5478.89)
+    SplineCoef.append(904.097)
+    SplineCoef.append(2.8)
+    SplineCoef.append(-115.387)
+    SplineCoef.append(888.81)
+    SplineCoef.append(-232.992)
+    SplineCoef.append(2.35)
+    SplineCoef.append(17587.6)
+    SplineCoef.append(-14177.6)
+    SplineCoef.append(2972.62)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Fullerene_RhoVsLogDspline_Var5(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
+    # units: rho [kg/m³]; diameter [nm]
+    # note: experimental mobility/mass data taken from: Moteki et al., Aerosol Sci. Technol., 44, 663-675, 2010;
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(13295.3)
+    SplineCoef.append(-8550.29)
+    SplineCoef.append(1423.82)
+    SplineCoef.append(2.8)
+    SplineCoef.append(-4979.54)
+    SplineCoef.append(4503.18)
+    SplineCoef.append(-907.154)
+    SplineCoef.append(2.45)
+    SplineCoef.append(8500)
+    SplineCoef.append(-6500.53)
+    SplineCoef.append(1338.5)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Fullerene_RhoVsLogDspline_Var8(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
+    # units: rho [kg/m³]; diameter [nm]
+    # note: experimental mobility/mass data taken from: PSI data for fullerene soot sample obtained from Shuka Schwarz, measured at ETH in Oct. 2010;
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(-2281.86)
+    SplineCoef.append(2687.72)
+    SplineCoef.append(-613.973)
+    SplineCoef.append(2.735)
+    SplineCoef.append(-1987.85)
+    SplineCoef.append(2472.72)
+    SplineCoef.append(-574.667)
+    SplineCoef.append(2.35)
+    SplineCoef.append(7319.9)
+    SplineCoef.append(-5448.76)
+    SplineCoef.append(1110.76)
+    SplineCoef.append(1.85)
+    SplineCoef.append(-520.822)
+    SplineCoef.append(3027.69)
+    SplineCoef.append(-1180.18)
+    SplineCoef.append(1.48)
+    SplineCoef.append(1421.03)
+    SplineCoef.append(403.564)
+    SplineCoef.append(-293.649)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Fullerene_RhoVsLogDspline_Var8_old(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
+    # units: rho [kg/m³]; diameter [nm]
+    # note: experimental mobility/mass data taken from: PSI data for fullerene soot sample obtained from Shuka Schwarz, measured at ETH in Oct. 2010;
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(6744.81)
+    SplineCoef.append(-4200.82)
+    SplineCoef.append(700)
+    SplineCoef.append(2.57)
+    SplineCoef.append(-6902.08)
+    SplineCoef.append(6419.32)
+    SplineCoef.append(-1366.18)
+    SplineCoef.append(2.4)
+    SplineCoef.append(7288.07)
+    SplineCoef.append(-5405.8)
+    SplineCoef.append(1097.39)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Aquadag_RhoVsLogDspline_Var25(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
+    # units: rho [kg/m³]; diameter [nm]
+    # note: experimental mobility/mass data from: Moteki et al., Aerosol Sci. Technol., 44, 663-675, 2010 (and they agree with pers. comm. from 04/09/2009)
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(10000)
+    SplineCoef.append(-6088.15)
+    SplineCoef.append(974.717)
+    SplineCoef.append(2.8)
+    SplineCoef.append(-3200)
+    SplineCoef.append(3340.42)
+    SplineCoef.append(-708.956)
+    SplineCoef.append(2.35)
+    SplineCoef.append(6490.76)
+    SplineCoef.append(-4907.03)
+    SplineCoef.append(1045.82)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Aquadag_RhoVsLogDspline_Var8(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
+    # units: rho [kg/m³]; diameter [nm]
+    # note: experimental mobility/mass data taken from: "grand average" of PSI's (APM at ETH in Oct. 2010) and DMT's (Subu; CPMA in 2011 by XXXX) data for classic and new aquadag batches;
+    if SplineCoef is None:
+        SplineCoef = []
+    SplineCoef.append(10092.9)
+    SplineCoef.append(-6171.83)
+    SplineCoef.append(970.172)
+    SplineCoef.append(2.75)
+    SplineCoef.append(-2627.45)
+    SplineCoef.append(3079.33)
+    SplineCoef.append(-711.856)
+    SplineCoef.append(2.25)
+    SplineCoef.append(6140.89)
+    SplineCoef.append(-4714.75)
+    SplineCoef.append(1020.16)
+    SplineCoef.append(1.65)
+    SplineCoef.append(1063.19)
+    SplineCoef.append(1440.04)
+    SplineCoef.append(-844.928)
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+def Aquadag_RhoVsLogDspline_Var8_old(
+    SplineCoef: list[float] | None = None,
+) -> tuple[list[float], int]:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function creates the spline fit coefficients for the relationship: density versus log(diam)
+    # units: rho [kg/m³]; diameter [nm]
+    # note: experimental mobility/mass data taken from: "grand average" of PSI's (APM at ETH in Oct. 2010) and DMT's (Subu; CPMA in 2011 by XXXX) data for classic and new aquadag batches;
+    if SplineCoef is None:
+        SplineCoef = [0] * 15  # Pre-initialize with 15 zeros for index assignment
+    SplineCoef[0] = 7900.57
+    SplineCoef[1] = -4964.2
+    SplineCoef[2] = 830.317
+    SplineCoef[3] = 2.6
+    SplineCoef[4] = -5371.15
+    SplineCoef[5] = 5244.82
+    SplineCoef[6] = -1132.96
+    SplineCoef[7] = 2.35
+    SplineCoef[8] = 5809.57
+    SplineCoef[9] = -4270.69
+    SplineCoef[10] = 891.62
+    SplineCoef[11] = 1.6
+    SplineCoef[12] = 904.355
+    SplineCoef[13] = 1860.83
+    SplineCoef[14] = -1024.48
+    return SplineCoef, SP2_calibCurveSplineCheck(SplineCoef)
+
+
+# Glassy Carbon functions
+def GlassyCarbonAlpha_Mass2Diam(mass: float | np.ndarray) -> float | np.ndarray:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    diam = 1e3 * (6 * mass / np.pi / 1420) ** (1 / 3)
+    return diam
+
+
+def GlassyCarbonAlpha_Diam2Mass(DiamMob: float | np.ndarray) -> float | np.ndarray:
+    # this function comes from the sp2 toolkit. Substantially it is a wrapper for the calib coeff
+    # This function calculates the mass of a glassy carbon particle of a given mobility diameter.
+    # exact description: Glassy carbon spherical powder, 0.4-12 micron, type 2, Alpha Aesar Art. No. 38008
+    # note: the bulk density is taken from the supplier's data sheet and an email by Ronald Becht (Technical service, Alfa Aesar)
+    kGlassyCarbonAlphaRho = 1420  # kg/m3
+    mass = kGlassyCarbonAlphaRho * np.pi / 6 * DiamMob**3 * 1e-9
+    return mass