From 6fa6fabf03e3c074382b72fb4f511e23ef1ddf84 Mon Sep 17 00:00:00 2001 From: Barbara Bertozzi Date: Wed, 6 Aug 2025 09:00:16 +0200 Subject: [PATCH] refactor: moved some functions from toolkit_legacy --- src/sp2xr/calibration.py | 449 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 449 insertions(+) diff --git a/src/sp2xr/calibration.py b/src/sp2xr/calibration.py index 7f22e0f..8bc9718 100644 --- a/src/sp2xr/calibration.py +++ b/src/sp2xr/calibration.py @@ -1,4 +1,5 @@ import numpy as np +from numpy.polynomial import Polynomial def polynomial(x, *coefficients): @@ -127,3 +128,451 @@ def calibrate_particle_data(df, config): ) return df + + +# %% Convertions mass to D , N to mass + + +def mass2meqDiam(mass_array, rho=1800): + """ + Calculate mass equivalent diameters for the input array of mass values. + inputs: mass_array values in fg, particle density in kg m-3 (rho) + outputs: array of mass equivalent diameters in nm + """ + Dp = np.array([1e9 * (6 * m * 1e-18 / rho / np.pi) ** (1 / 3) for m in mass_array]) + return Dp + + +def dNdlogDp_to_dMdlogDp(N, Dp, rho=1800): + """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] + return the corresponding the mass size distribution [dM/dlogDp, ug/m3] + + Optional density input given in units of kg/m3 + """ + # M = N*((np.pi*(Dp*1e-9)**3)/6)*rho*1e18 + M = N * np.pi / 6 * rho * (Dp**3) * 1e-12 + return M + + +def BC_mass_to_diam(mass, rho_eff=1800, BC_type=""): + """ + Thi function calculates the effective density and mobility (volume equivalent?) diameter + for BC calibration materials applying different parameterizations. + + Parameters + ---------- + mass : float or array of floats + BC mass in fg. + rho_eff : float, optional + Effective density [kg/m3] to use in case BC_type='constant_mobility_density'. The default is 1800. + BC_type : string, optional + Identification string for the effective density parameterization to use. The default is ''. + + Raises + ------ + ValueError + If the string provided in BC_type is not inlcuded in the list below an error is raised. + + Returns + ------- + density : float of array of floats + Effective density [kg/m3] corresponding to the input mass(es). + diam : float or array of floats + Diameter [nm] corresponding to the mass and effective density. + + """ + if BC_type == "constant_effective_density": + density = np.zeros_like(mass) + rho_eff + diam = 1e3 * ((6 * mass) / (np.pi * rho_eff)) ** (1 / 3) + elif BC_type == "Aquadag_RCAST_2009": + SplineCoefM2D, nSplineSegmM2D = Aquadag_RhoVsLogMspline_Var25(SplineCoef=[]) + density = SP2_calibCurveSpline(mass, SplineCoefM2D) + elif BC_type == "Aquadag_Moteki_AST2010": + SplineCoefM2D, nSplineSegmM2D = Aquadag_RhoVsLogMspline_Var25( + SplineCoef=[] + ) # this line only loads the coefficients from Moteki 2010 + density = SP2_calibCurveSpline(mass, SplineCoefM2D) + elif BC_type == "Aquadag_PSI_2010": + SplineCoefM2D, nSplineSegmM2D = Aquadag_RhoVsLogMspline_Var8(SplineCoef=[]) + density = SP2_calibCurveSpline(mass, SplineCoefM2D) + elif BC_type == "FS_RCAST_2009": + SplineCoefM2D, nSplineSegmM2D = Fullerene_RhoVsLogMspline_Var2(SplineCoef=[]) + density = SP2_calibCurveSpline(mass, SplineCoefM2D) + elif BC_type == "FS_Moteki_2010": + SplineCoefM2D, nSplineSegmM2D = Fullerene_RhoVsLogMspline_Var5(SplineCoef=[]) + density = SP2_calibCurveSpline(mass, SplineCoefM2D) + elif BC_type == "FS_PSI_2010": + SplineCoefM2D, nSplineSegmM2D = Fullerene_RhoVsLogMspline_Var8(SplineCoef=[]) + density = SP2_calibCurveSpline(mass, SplineCoefM2D) + elif BC_type == "FS_PSI_2010_old": + SplineCoefM2D, nSplineSegmM2D = Fullerene_RhoVsLogMspline_Var8_old( + SplineCoef=[] + ) + density = SP2_calibCurveSpline(mass, SplineCoefM2D) + elif BC_type == "Glassy_Carbon": + diam = GlassyCarbonAlpha_Mass2Diam(mass) + else: + raise ValueError("Calibration material not defined") + + diam = 1e3 * (6 * mass / np.pi / density) ** (1 / 3) # 1e3: unit conversion + + return density, diam + + +def BC_diam_to_mass(diam, rho_eff=1800, BC_type=""): + + if BC_type == "": # pers comm Kondo 2009, spline fit + SplineCoefD2M, nSplineSegmD2M = Fullerene_RhoVsLogDspline_Var2(SplineCoef=[]) + density = SP2_calibCurveSpline(diam, SplineCoefD2M) + elif BC_type == "FS_Moteki_2010": # Moteki et al., AST, 2010 + SplineCoefD2M, nSplineSegmD2M = Fullerene_RhoVsLogDspline_Var5(SplineCoef=[]) + density = SP2_calibCurveSpline(diam, SplineCoefD2M) + elif BC_type == "FS_PSI_2010": + SplineCoefD2M, nSplineSegmD2M = Fullerene_RhoVsLogDspline_Var8(SplineCoef=[]) + density = SP2_calibCurveSpline(diam, SplineCoefD2M) + elif BC_type == "FS_PSI_2010_old": + SplineCoefD2M, nSplineSegmD2M = Fullerene_RhoVsLogDspline_Var8_old( + SplineCoef=[] + ) + density = SP2_calibCurveSpline(diam, SplineCoefD2M) + + else: + raise ValueError("Calibration material not defined") + + 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)