diff --git a/carculator_truck/__init__.py b/carculator_truck/__init__.py
index 552d650..6066b75 100644
--- a/carculator_truck/__init__.py
+++ b/carculator_truck/__init__.py
@@ -22,7 +22,7 @@
 )
 
 # library version
-__version__ = (0, 2, 6)
+__version__ = (0, 2, 7)
 
 from pathlib import Path
 
diff --git a/carculator_truck/data/A_matrix.csv b/carculator_truck/data/A_matrix.csv
index b65ec92..bb7fc95 100644
--- a/carculator_truck/data/A_matrix.csv
+++ b/carculator_truck/data/A_matrix.csv
@@ -1743,7 +1743,7 @@
 534;939;6.89E-07
 535;535;1
 536;536;1
-536;915;4
+536;915;3.57
 537;537;1
 538;39;-42840
 538;347;-0.225
@@ -2534,7 +2534,7 @@
 772;819;-1
 772;850;-0.059579832
 772;886;-0.01304
-772;915;-3
+772;915;-0.8
 772;922;-0.1571
 772;957;-0.15
 773;80;-0.00030042
@@ -3246,3 +3246,4 @@
 985;985;1
 986;408;-0.6
 986;986;1
+536;915;-1.77
diff --git a/carculator_truck/data/extra_parameters.json b/carculator_truck/data/extra_parameters.json
index 738b2cf..3683533 100644
--- a/carculator_truck/data/extra_parameters.json
+++ b/carculator_truck/data/extra_parameters.json
@@ -259,5 +259,8 @@
   "Mercury direct emissions, rural",
   "Cadmium direct emissions, rural",
   "battery lifetime kilometers",
-  "capacity utilization"
+  "capacity utilization",
+  "share recuperated energy",
+  "road dust emissions",
+  "transmission efficiency"
 ]
diff --git a/carculator_truck/energy_consumption.py b/carculator_truck/energy_consumption.py
index 31002de..71575f0 100644
--- a/carculator_truck/energy_consumption.py
+++ b/carculator_truck/energy_consumption.py
@@ -155,13 +155,14 @@ def motive_energy_per_km(
         gradient_resistance = (driving_mass * 9.81).T.values * np.sin(self.gradient)
         # Inertia: driving mass * acceleration
         inertia = self.acceleration * driving_mass.values.T
-        # Braking loss: when inertia is negative
-        braking_loss = np.where(inertia < 0, inertia * -1, 0)
 
         total_resistance = (
             rolling_resistance + air_resistance + gradient_resistance + inertia
         )
 
+        # Braking loss: when inertia is negative
+        braking_loss = np.where(total_resistance < 0, total_resistance * -1, 0)
+
         if not debug_mode:
 
             # Power required: total resistance * velocity
diff --git a/carculator_truck/inventory.py b/carculator_truck/inventory.py
index 195a3fa..c2651b3 100644
--- a/carculator_truck/inventory.py
+++ b/carculator_truck/inventory.py
@@ -1670,19 +1670,12 @@ def add_additional_activities(self):
 
                     if pt == "BEV":
 
-                        name = (
-                            "transport, freight, lorry, "
-                            + pt
-                            + ", "
-                            + self.background_configuration["energy storage"][
+                        chemistry = self.background_configuration["energy storage"][
                                 "electric"
                             ]["BEV"]
-                            + " battery, "
-                            + s
-                            + " gross weight, "
-                            + str(y)
-                            + ", "
-                            + self.cycle.lower()
+
+                        name = (
+                            f"transport, freight, lorry, {pt}, {chemistry} battery, {s} gross weight, {y}, {self.cycle.lower()}"
                         )
 
                         if self.scope["fu"]["unit"] == "tkm":
@@ -1701,6 +1694,10 @@ def add_additional_activities(self):
 
                     elif pt == "FCEV":
 
+                        name = (
+                            f"transport, freight, lorry, {pt}, {s} gross weight, {y}, {self.cycle.lower()}"
+                        )
+
                         name = (
                             "transport, freight, lorry, "
                             + pt
@@ -1729,16 +1726,7 @@ def add_additional_activities(self):
                     else:
 
                         name = (
-                            "transport, freight, lorry, "
-                            + pt
-                            + ", "
-                            + s
-                            + " gross weight, "
-                            + str(y)
-                            + ", "
-                            + euro_class
-                            + ", "
-                            + self.cycle.lower()
+                            f"transport, freight, lorry, {pt}, {s} gross weight, {y}, {euro_class}, {self.cycle.lower()}"
                         )
 
                         if self.scope["fu"]["unit"] == "tkm":
@@ -1751,7 +1739,7 @@ def add_additional_activities(self):
                                 name,
                                 self.country,
                                 unit,
-                                "transport, freight, lorry, " + euro_class,
+                                f"transport, freight, lorry, {euro_class}",
                             )
                         ] = maximum
 
@@ -5088,9 +5076,6 @@ def set_inputs_in_A_matrix(self, array):
         )
 
         # Brake wear emissions
-        # BEVs only emit 20% of what a combustion engine vehicle emit according to
-        # https://link.springer.com/article/10.1007/s11367-014-0792-4
-
         self.A[
             :,
             self.inputs[
@@ -5106,27 +5091,6 @@ def set_inputs_in_A_matrix(self, array):
             array[self.array_inputs["total cargo mass"], :] / 1000
         )
 
-        ind_A = [
-            i
-            for i in self.car_indices
-            if any(
-                x in self.rev_inputs[i][0] for x in ["BEV", "FCEV", "HEV-d", "PHEV-d"]
-            )
-        ]
-
-        self.A[
-            :,
-            self.inputs[
-                (
-                    "treatment of brake wear emissions, lorry",
-                    "RER",
-                    "kilogram",
-                    "brake wear emissions, lorry",
-                )
-            ],
-            ind_A,
-        ] *= 0.2
-
         # Infrastructure: 5.37e-4 per gross tkm
         self.A[
             :,
@@ -6657,10 +6621,6 @@ def set_inputs_in_A_matrix_for_export(self, array):
         )
 
         # Brake wear emissions
-        # BEVs and other hybrid vehicles only emit 20%
-        # of what a combustion engine vehicle emit according to
-        # https://link.springer.com/article/10.1007/s11367-014-0792-4
-
         self.A[
             :,
             self.inputs[
@@ -6676,27 +6636,6 @@ def set_inputs_in_A_matrix_for_export(self, array):
             array[self.array_inputs["total cargo mass"], :] / 1000
         )
 
-        ind_A = [
-            i
-            for i in self.car_indices
-            if any(
-                x in self.rev_inputs[i][0] for x in ["BEV", "FCEV", "HEV-d", "PHEV-d"]
-            )
-        ]
-
-        self.A[
-            :,
-            self.inputs[
-                (
-                    "treatment of brake wear emissions, lorry",
-                    "RER",
-                    "kilogram",
-                    "brake wear emissions, lorry",
-                )
-            ],
-            ind_A,
-        ] *= 0.2
-
         # Infrastructure: 5.37e-4 per gross tkm
         self.A[
             :,
diff --git a/carculator_truck/model.py b/carculator_truck/model.py
index cea0467..b57da95 100644
--- a/carculator_truck/model.py
+++ b/carculator_truck/model.py
@@ -12,12 +12,7 @@
 
 DEFAULT_MAPPINGS = {
     "electric": {"BEV", "PHEV-e"},
-    "combustion": {
-        "HEV-d",
-        "ICEV-g",
-        "ICEV-d",
-        "PHEV-c-d",
-    },
+    "combustion": {"HEV-d", "ICEV-g", "ICEV-d", "PHEV-c-d",},
     "combustion_wo_cng": {"HEV-d", "ICEV-d", "PHEV-c-d"},
     "pure_combustion": {"ICEV-g", "ICEV-d"},
     "petrol": {"PHEV-c-p"},
@@ -134,9 +129,7 @@ def __init__(
                             if s in self.array.coords["size"].values
                         ],
                     )
-                ] *= (
-                    1 - 0.67
-                )
+                ] *= (1 - 0.67)
 
             if any(
                 s in self.array.coords["size"].values for s in ["32t", "40t", "60t"]
@@ -150,9 +143,7 @@ def __init__(
                             if s in self.array.coords["size"].values
                         ],
                     )
-                ] *= (
-                    1 - 0.33
-                )
+                ] *= (1 - 0.33)
 
         if self.cycle == "Urban delivery":
             self.array.loc[dict(parameter="kilometers per year")] *= 1 - 0.39
@@ -170,9 +161,7 @@ def __init__(
                             if s in self.array.coords["size"].values
                         ],
                     )
-                ] *= (
-                    1 - 0.67
-                )
+                ] *= (1 - 0.67)
 
             if any(
                 s in self.array.coords["size"].values for s in ["32t", "40t", "60t"]
@@ -186,9 +175,7 @@ def __init__(
                             if s in self.array.coords["size"].values
                         ],
                     )
-                ] *= (
-                    1 - 0.33
-                )
+                ] *= (1 - 0.33)
 
     def __call__(self, key):
         """
@@ -278,6 +265,7 @@ def set_all(self):
             self.set_auxiliaries()
             self.set_fuel_cell_parameters()
             self.calculate_ttw_energy()
+            self.set_share_recuperated_energy()
             self.set_battery_fuel_cell_replacements()
             self.set_energy_stored_properties()
 
@@ -300,11 +288,14 @@ def set_all(self):
         self.set_electric_utility_factor()
         self.set_electricity_consumption()
         self.set_costs()
+        self.create_PHEV()
+        self.drop_hybrid()
         self.set_hot_emissions()
         self.set_particulates_emission()
         self.set_noise_emissions()
-        self.create_PHEV()
-        self.drop_hybrid()
+
+        self.array.values = np.clip(self.array.values, 0, None)
+
 
         print("")
         print("'-' BEV with driving mass superior to the permissible gross weight.")
@@ -455,16 +446,13 @@ def adjust_combustion_power_share(self):
                 )
 
                 if arr.sum() > 0:
-                    new_shares = (
-                        self.array.loc[
-                            dict(
-                                powertrain=l_pwt,
-                                parameter="combustion power share",
-                                year=actual_years,
-                            )
-                        ]
-                        - (arr * 0.02)
-                    )
+                    new_shares = self.array.loc[
+                        dict(
+                            powertrain=l_pwt,
+                            parameter="combustion power share",
+                            year=actual_years,
+                        )
+                    ] - (arr * 0.02)
                     self.array.loc[
                         dict(
                             powertrain=l_pwt,
@@ -552,11 +540,7 @@ def adjust_cost(self):
         # Correction of combustion powertrain cost for ICEV-g
         if "ICEV-g" in self.array.powertrain.values:
             self.array.loc[
-                :,
-                ["ICEV-g"],
-                "combustion powertrain cost per kW",
-                :,
-                :,
+                :, ["ICEV-g"], "combustion powertrain cost per kW", :, :,
             ] = np.reshape(
                 (5.92e160 * np.exp(-0.1819 * self.array.year.values) + 26.76)
                 * cost_factor,
@@ -572,15 +556,7 @@ def drop_hybrid(self):
         l_pwt = [
             p
             for p in self.array.powertrain.values
-            if p
-            in [
-                "ICEV-d",
-                "ICEV-g",
-                "PHEV-d",
-                "FCEV",
-                "BEV",
-                "HEV-d",
-            ]
+            if p in ["ICEV-d", "ICEV-g", "PHEV-d", "FCEV", "BEV", "HEV-d",]
         ]
         self.array = self.array.sel(powertrain=l_pwt)
 
@@ -608,7 +584,7 @@ def calculate_ttw_energy(self):
                 (
                     len(self.array.coords["size"]),
                     len(self.array.coords["powertrain"]),
-                    8,
+                    9,
                     len(self.array.coords["year"]),
                     len(self.array.coords["value"]),
                     self.ecm.cycle.shape[0],
@@ -626,6 +602,7 @@ def calculate_ttw_energy(self):
                     "transmission efficiency",
                     "engine efficiency",
                     "power load",
+                    "negative motive energy",
                 ],
                 self.array.coords["year"],
                 self.array.coords["value"],
@@ -651,6 +628,10 @@ def calculate_ttw_energy(self):
             motive_power.T, 0, self["power"].values[..., None]
         )
 
+        self.energy.loc[dict(parameter="negative motive energy")] = (
+            np.clip(motive_power.T, None, 0) * -1
+        )
+
         self.energy.loc[dict(parameter="recuperated energy at wheels")] = (
             np.clip(recuperated_power.T, 0, self["electric power"].values[..., None])
             * -1
@@ -669,10 +650,11 @@ def calculate_ttw_energy(self):
         if len(l_pwt) > 0:
 
             self.energy.loc[dict(parameter="power load", powertrain=l_pwt)] = (
-                motive_power.T[:, idx, ...] + recuperated_power.T[:, idx, ...]
-            ) / self.array.sel(parameter="electric power", powertrain=l_pwt).values[
-                ..., None
-            ]
+                (motive_power.T[:, idx, ...] + recuperated_power.T[:, idx, ...])
+                / self.array.sel(parameter="electric power", powertrain=l_pwt).values[
+                    ..., None
+                ]
+            )
 
         self.energy.loc[dict(parameter="power load")] = np.clip(
             self.energy.loc[dict(parameter="power load")], 0, 1
@@ -697,10 +679,7 @@ def calculate_ttw_energy(self):
         if len(l_pwt) > 0:
 
             self.energy.loc[
-                dict(
-                    parameter="engine efficiency",
-                    powertrain=l_pwt,
-                )
+                dict(parameter="engine efficiency", powertrain=l_pwt,)
             ] = np.clip(
                 np.interp(
                     self.energy.loc[dict(parameter="power load", powertrain=l_pwt)],
@@ -728,35 +707,25 @@ def calculate_ttw_energy(self):
         if "ICEV-g" in self.array.powertrain.values:
             self.energy.loc[
                 dict(parameter="engine efficiency", powertrain="ICEV-g")
-            ] *= 1 - self.array.sel(
-                parameter="CNG engine efficiency correction factor", powertrain="ICEV-g"
+            ] *= (
+                1
+                - self.array.sel(
+                    parameter="CNG engine efficiency correction factor",
+                    powertrain="ICEV-g",
+                )
             )
 
-        # Correction for electric motors
-
-        l_pwt = [p for p in self.array.powertrain.values if p in ["BEV", "PHEV-e"]]
-
+        l_pwt = [
+            p
+            for p in self.array.powertrain.values
+            if p in ["BEV", "FCEV", "PHEV-e"]
+        ]
         if len(l_pwt) > 0:
-
-            self.energy.loc[dict(parameter="engine efficiency", powertrain=l_pwt)] = (
-                self.array.sel(parameter="engine efficiency", powertrain=l_pwt)
-                * self.array.sel(
-                    parameter="battery discharge efficiency", powertrain=l_pwt
-                )
-            ).expand_dims(dim="x", axis=-1)
-
-        # Correction for fuel cell stack efficiency
-        if "FCEV" in self.array.powertrain.values:
             self.energy.loc[
-                dict(parameter="engine efficiency", powertrain=["FCEV"])
-            ] = (
-                self.array.sel(
-                    parameter="fuel cell system efficiency", powertrain=["FCEV"]
-                )
-                * self.array.sel(parameter="engine efficiency", powertrain=["FCEV"])
-            ).expand_dims(
-                dim="x", axis=-1
-            )
+                dict(parameter="engine efficiency", powertrain=l_pwt)
+            ] = self.array.loc[dict(
+                parameter="engine efficiency", powertrain=l_pwt
+            )].values[..., None]
 
         self.energy.loc[dict(parameter="motive energy")] = (
             self.energy.loc[dict(parameter="motive energy at wheels")]
@@ -798,17 +767,37 @@ def calculate_ttw_energy(self):
             / distance
         ).T
 
-        self["TtW efficiency"] = (
-            (
+        self["TtW efficiency"] = np.nanmean(
+            np.where(
+                self.energy.loc[dict(parameter="power load")] == 0,
+                np.nan,
                 (
-                    self.energy.sel(parameter="motive energy at wheels").sum(
-                        dim="second"
-                    )
-                    + self.energy.sel(parameter="auxiliary energy").sum(dim="second")
-                ).T
-                / distance
-            ).T
-        ) / self["TtW energy"]
+                    self.energy.loc[dict(parameter="transmission efficiency")]
+                    * self.energy.loc[dict(parameter="engine efficiency")]
+                ),
+            ),
+            axis=-1,
+        )
+
+        # Correction for electric motors
+        l_pwt = [p for p in self.array.powertrain.values if p in ["BEV", "PHEV-e"]]
+
+        if len(l_pwt) > 0:
+            self.array.loc[
+                dict(parameter="TtW efficiency", powertrain=l_pwt)
+            ] *= self.array.sel(
+                parameter="battery discharge efficiency", powertrain=l_pwt
+            )
+
+        # Correction for fuel cell stack efficiency
+        if "FCEV" in self.array.powertrain.values:
+            self.array.loc[
+                dict(parameter="TtW efficiency", powertrain=["FCEV"])
+            ] *= self.array.sel(
+                parameter="fuel cell system efficiency", powertrain=["FCEV"]
+            ) * self.array.sel(
+                parameter="battery discharge efficiency", powertrain=["FCEV"]
+            )
 
         self["auxiliary energy"] = (
             self.energy.sel(parameter="auxiliary energy").sum(dim="second").T / distance
@@ -831,6 +820,15 @@ def calculate_ttw_energy(self):
             axis=-1,
         )
 
+        self.array.loc[dict(parameter="transmission efficiency")] = np.nanmean(
+            np.where(
+                self.energy.loc[dict(parameter="power load")] == 0,
+                np.nan,
+                self.energy.loc[dict(parameter="transmission efficiency")],
+            ),
+            axis=-1,
+        )
+
     def set_fuel_cell_parameters(self):
         """
         Specific setup for fuel cells, which are mild hybrids.
@@ -1064,6 +1062,27 @@ def set_component_masses(self):
             + self["inverter fix mass"]
         )
 
+    def set_share_recuperated_energy(self):
+        """ Calculate the share of recuperated energy, over the total negative motive energy"""
+
+        self["share recuperated energy"] = (
+            (
+                self.energy.loc[dict(parameter="recuperated energy at wheels")].sum(
+                    dim="second"
+                )
+                * -1
+            )
+            * self["engine efficiency"]
+            * self["transmission efficiency"]
+        ) / self.energy.loc[dict(parameter="negative motive energy")].sum(dim="second")
+
+        if "PHEV-d" in self.array.powertrain.values:
+            self.array.loc[
+                dict(powertrain="PHEV-c-d", parameter="share recuperated energy")
+            ] = self.array.loc[
+                dict(powertrain="PHEV-e", parameter="share recuperated energy")
+            ]
+
     def set_electric_utility_factor(self):
         """
         The electric utility factor
@@ -1107,12 +1126,83 @@ def create_PHEV(self):
                 )
             )
 
+            #self.array.loc[{"powertrain": "PHEV-d", "parameter":[
+            #    "battery cell mass", "battery BoP mass", "energy battery mass"
+            #]}] = self.array.loc[{"powertrain": "PHEV-e", "parameter":[
+            #    "battery cell mass", "battery BoP mass", "energy battery mass"
+            #]}]
+
+
             self.array.loc[
                 {"powertrain": "PHEV-d", "parameter": "electric utility factor"}
             ] = self.array.loc[
                 {"powertrain": "PHEV-e", "parameter": "electric utility factor"}
             ]
 
+            self.energy.loc[
+                dict(
+                    parameter=[
+                        "auxiliary energy",
+                        "motive energy",
+                        "motive energy at wheels",
+                        "recuperated energy",
+                        "recuperated energy at wheels",
+                        "transmission efficiency",
+                        "engine efficiency",
+                        "power load",
+                        "negative motive energy",
+                    ],
+                    powertrain=["PHEV-d"],
+                )
+            ] = (
+                self.array.loc[
+                    dict(parameter="electric utility factor", powertrain=["PHEV-e"])
+                ]
+                * self.energy.loc[
+                    dict(
+                        parameter=[
+                            "auxiliary energy",
+                            "motive energy",
+                            "motive energy at wheels",
+                            "recuperated energy",
+                            "recuperated energy at wheels",
+                            "transmission efficiency",
+                            "engine efficiency",
+                            "power load",
+                            "negative motive energy",
+                        ],
+                        powertrain=["PHEV-e"],
+                    )
+                ]
+            ).values.transpose(
+                0, 1, 4, 2, 3, 5
+            ) + (
+                (
+                    1
+                    - self.array.loc[
+                        dict(parameter="electric utility factor", powertrain="PHEV-e")
+                    ]
+                )
+                * self.energy.loc[
+                    dict(
+                        parameter=[
+                            "auxiliary energy",
+                            "motive energy",
+                            "motive energy at wheels",
+                            "recuperated energy",
+                            "recuperated energy at wheels",
+                            "transmission efficiency",
+                            "engine efficiency",
+                            "power load",
+                            "negative motive energy",
+                        ],
+                        powertrain=["PHEV-c-d"],
+                    )
+                ]
+            ).values.transpose(
+                0, 3, 4, 1, 2, 5
+            )
+
     def set_energy_stored_properties(self):
         """
         First, fuel mass is defined. It is dependent on the range required.
@@ -1402,17 +1492,27 @@ def set_particulates_emission(self):
         Calculate the emission of particulates according to
         https://www.eea.europa.eu/ds_resolveuid/6USNA27I4D
 
+        and further disaggregated in:
+        https://doi.org/10.1016/j.atmosenv.2020.117886
+
         for:
 
         - brake wear
         - tire wear
-        - road wera
+        - road wear
+        - re-suspended road dust
+
+        by considering:
+
+        - vehicle mass
+        - driving situation (urban, rural, motorway)
 
-        Tier-2 method considers:
+        into the following fractions:
 
-        - number of axles
-        - capacity utilization rate
-        - speed
+        - PM 2.5
+        - PM 10
+
+        Emissions are subdivided in compartments: urban, suburban and rural.
 
         """
 
@@ -1420,25 +1520,20 @@ def set_particulates_emission(self):
             "tire wear emissions",
             "brake wear emissions",
             "road wear emissions",
+            "road dust emissions",
         ]
 
         pem = ParticulatesEmissionsModel(
             cycle_name=self.ecm.cycle_name,
             cycle=self.ecm.cycle,
-            number_axles=self["number of axles"],
-            load_factor=self["capacity utilization"],
+            mass=self["driving mass"],
         )
 
         res = pem.get_abrasion_emissions()
-        self[list_param] = res.transpose(4, 3, 0, 2, 1)
+        self[list_param] = res
 
-        # Superseed values for 3.5t and 7.5t trucks
-        for size in [
-            s for s in ["3.5t", "7.5t"] if s in self.array.coords["size"].values
-        ]:
-            self.array.loc[dict(size=[size], parameter=list_param)] = (
-                np.array([0.0169, 0.0117, 0.015])[None, :, None, None] / 1000
-            )
+        # brake emissions are discounted by the use of regenerative braking
+        self["brake wear emissions"] *= 1 - self["share recuperated energy"]
 
     def set_hot_emissions(self):
         """
@@ -1588,18 +1683,12 @@ def set_hot_emissions(self):
 
         if len(l_pwt) > 0:
             self.array.loc[
-                dict(
-                    powertrain=l_pwt,
-                    parameter=list_direct_emissions,
-                )
+                dict(powertrain=l_pwt, parameter=list_direct_emissions,)
             ] = hem.get_emissions_per_powertrain(
                 powertrain_type="diesel",
                 euro_classes=l_y,
                 lifetime_km=self.array.loc[
-                    dict(
-                        powertrain=l_pwt,
-                        parameter="lifetime kilometers",
-                    )
+                    dict(powertrain=l_pwt, parameter="lifetime kilometers",)
                 ],
                 energy_consumption=self.energy.sel(
                     powertrain=l_pwt,
@@ -1620,10 +1709,7 @@ def set_hot_emissions(self):
                     powertrain_type="cng",
                     euro_classes=l_y,
                     lifetime_km=self.array.loc[
-                        dict(
-                            powertrain="ICEV-g",
-                            parameter="lifetime kilometers",
-                        )
+                        dict(powertrain="ICEV-g", parameter="lifetime kilometers",)
                     ],
                     energy_consumption=self.energy.sel(
                         powertrain=["ICEV-g"],
@@ -1858,10 +1944,7 @@ def get_share_biofuel(self):
 
         share_biofuel = (
             self.bs.biofuel.sel(
-                region=region,
-                value=0,
-                fuel_type="Biomass fuel",
-                scenario=scenario,
+                region=region, value=0, fuel_type="Biomass fuel", scenario=scenario,
             )
             .interp(
                 year=self.array.coords["year"].values,
diff --git a/carculator_truck/particulates_emissions.py b/carculator_truck/particulates_emissions.py
index de6bcde..5460e02 100644
--- a/carculator_truck/particulates_emissions.py
+++ b/carculator_truck/particulates_emissions.py
@@ -3,27 +3,28 @@
 
 class ParticulatesEmissionsModel:
     """
-    Calculate particulates emissions based on the method described in
+    Calculate particulates emissions based on the method described in:
     https://www.eea.europa.eu/ds_resolveuid/6USNA27I4D
 
+    and further disaggregated in:
+    https://doi.org/10.1016/j.atmosenv.2020.117886
+
     Include emission from:
 
     - brake wear
     - tire wear
     - road wear
+    - re-suspended road dust
 
     by considering:
 
-    - load factor of the vehicle
-    - number of axles
-    - speed levels
+    - vehicle mass
+    - driving situation (urban, rural, motorway)
 
     into the following fractions:
 
-    - PM 10
     - PM 2.5
-    - PM 1
-    - PM 0.1
+    - PM 10
 
     Emissions are subdivided in compartments: urban, suburban and rural.
 
@@ -34,16 +35,16 @@ class ParticulatesEmissionsModel:
 
     """
 
-    def __init__(self, cycle_name, cycle, number_axles, load_factor):
+    def __init__(self, cycle_name, cycle, mass):
 
         self.cycle_name = cycle_name
         self.cycle = np.resize(
             cycle,
             (
                 cycle.shape[0],
-                number_axles.shape[-1],
-                number_axles.shape[2],
-                number_axles.shape[1],
+                mass.shape[-1],
+                mass.shape[2],
+                mass.shape[1],
                 cycle.shape[-1],
             ),
         )
@@ -65,78 +66,135 @@ def __init__(self, cycle_name, cycle, number_axles, load_factor):
 
         self.velocity = self.cycle / 3600  # m/s per second
 
-        self.number_axles = number_axles.values
-        self.load_factor = load_factor.values
+        self.mass = mass.values / 1000  # in tons
 
     def get_abrasion_emissions(self):
 
-        tire_wear = (
-            self.get_tire_wear_emissions(self.number_axles, self.load_factor) / 1000
-        )
+        (
+            tire_pm10_urban,
+            tire_pm10_rural,
+            tire_pm10_motorway,
+            tire_pm25_urban,
+            tire_pm25_rural,
+            tire_pm25_motorway,
+        ) = self.get_tire_wear_emissions()
+
+        (
+            brake_pm10_urban,
+            brake_pm10_rural,
+            brake_pm10_motorway,
+            brake_pm25_urban,
+            brake_pm25_rural,
+            brake_pm25_motorway,
+        ) = self.get_brake_wear_emissions()
+
+        road_pm10, road_pm25 = self.get_road_wear_emissions()
 
-        brake_wear = self.get_brake_wear_emissions(self.load_factor) / 1000
+        dust_pm10, dust_pm25 = self.get_resuspended_road_dust()
 
-        road_wear = self.get_road_wear_emissions() / 1000
+        if self.cycle_name == "Urban delivery":
+            tire_wear = tire_pm10_urban + tire_pm25_urban
+            brake_wear = brake_pm10_urban + brake_pm25_urban
+
+        elif self.cycle_name == "Regional delivery":
+            tire_wear = (tire_pm10_urban + tire_pm25_urban) * 0.16
+            tire_wear += (tire_pm10_rural + tire_pm25_rural) * 0.32
+            tire_wear += (tire_pm10_motorway + tire_pm25_motorway) * 0.52
+
+            brake_wear = (brake_pm10_urban + brake_pm25_urban) * 0.16
+            brake_wear += (brake_pm10_rural + brake_pm25_rural) * 0.32
+            brake_wear += (brake_pm10_motorway + brake_pm25_motorway) * 0.52
+
+        else:
+            tire_wear = tire_pm10_motorway + tire_pm25_motorway
+            brake_wear = brake_pm10_motorway + brake_pm25_motorway
+
+        road_wear = road_pm10 + road_pm25
+
+        road_dust = dust_pm10 + dust_pm25
 
         res = np.vstack(
             (
-                tire_wear.sum(axis=0)[None, ...],
-                brake_wear.sum(axis=0)[None, ...],
-                road_wear.sum(axis=0)[None, ...],
+                tire_wear.sum(axis=2)[None, ...],
+                brake_wear.sum(axis=2)[None, ...],
+                road_wear.sum(axis=2)[None, ...],
+                road_dust.sum(axis=2)[None, ...],
             )
         )
 
-        distance = self.cycle.sum(axis=0) / 3600
-
-        return res / distance  # total emission per duty cycle to total emissions per km
+        return res.transpose(1, 2, 0, 3)[..., None]
 
-    def get_tire_wear_emissions(self, number_axles, load_factor):
+    def get_tire_wear_emissions(self):
         """
         Returns tire wear emissions.
 
-        :param number_axles: number of axles. Int.
-        :param load_factor: load factor. Float.
         :return:
         """
-        PM_total = number_axles / 2 * (1.41 + (1.38 * load_factor)) * 0.0107
 
-        PM_total = PM_total.T * self.velocity
-        PM_total[self.cycle <= 40] *= 1.39
-        PM_total[(self.cycle > 40) & (self.cycle < 90)] *= (
-            -0.00974 * self.cycle[(self.cycle > 40) & (self.cycle < 90)] + 1.78
-        )
-        PM_total[self.cycle > 90] *= 0.902
+        # converted to kg per vkm
+        pm10_urban = 5.8 * np.power(self.mass, 1 / 2.3) / 1e6
+        pm25_urban = 8.2 * np.power(self.mass, 1 / 2.3) / 1e6
+
+        pm10_rural = 4.5 * np.power(self.mass, 1 / 2.3) / 1e6
+        pm25_rural = 6.4 * np.power(self.mass, 1 / 2.3) / 1e6
+
+        pm10_motorway = 3.8 * np.power(self.mass, 1 / 2.3) / 1e6
+        pm25_motorway = 5.5 * np.power(self.mass, 1 / 2.3) / 1e6
 
-        return PM_total
+        return (
+            pm10_urban,
+            pm10_rural,
+            pm10_motorway,
+            pm25_urban,
+            pm25_rural,
+            pm25_motorway,
+        )
 
-    def get_brake_wear_emissions(self, load_factor):
+    def get_brake_wear_emissions(self):
         """
         Returns brake wear emissions.
 
-        :param number_axles: number of axles. Int.
-        :param load_factor: load factor. Float.
         :return:
         """
-        PM_total = 3.13 * (1 + 0.79 * load_factor) * 0.0075
-
-        PM_total = PM_total.T * self.velocity
-        PM_total[self.cycle <= 40] *= 1.67
-        PM_total[(self.cycle > 40) & (self.cycle < 90)] *= (
-            -0.027 * self.cycle[(self.cycle > 40) & (self.cycle < 90)] + 2.75
+        # converted to kg per vkm
+        pm10_urban = 4.2 * np.power(self.mass, 1 / 1.9) / 1e6
+        pm25_urban = 11 * np.power(self.mass, 1 / 1.9) / 1e6
+
+        pm10_rural = 1.8 * np.power(self.mass, 1 / 1.5) / 1e6
+        pm25_rural = 4.5 * np.power(self.mass, 1 / 1.5) / 1e6
+
+        pm10_motorway = 0.4 * np.power(self.mass, 1 / 1.3) / 1e6
+        pm25_motorway = 1.0 * np.power(self.mass, 1 / 1.3) / 1e6
+
+        return (
+            pm10_urban,
+            pm10_rural,
+            pm10_motorway,
+            pm25_urban,
+            pm25_rural,
+            pm25_motorway,
         )
-        PM_total[self.cycle > 90] *= 0.185
-
-        return PM_total
 
     def get_road_wear_emissions(self):
         """
         Returns road wear emissions.
 
-        :param number_axles: number of axles. Int.
-        :param load_factor: load factor. Float.
         :return:
         """
-        PM_total = np.full_like(self.velocity.T, 0.076)
-        PM_total = PM_total.T * self.velocity
+        # converted to kg per vkm
+        pm10 = 2.8 * np.power(self.mass, 1 / 1.5) / 1e6
+        pm25 = 5.1 * np.power(self.mass, 1 / 1.5) / 1e6
+
+        return (pm10, pm25)
+
+    def get_resuspended_road_dust(self):
+        """
+        Returns re-suspended road dust emissions.
+
+        :return:
+        """
+        # converted to kg per vkm
+        pm10 = 2 * np.power(self.mass, 1 / 1.1) / 1e6
+        pm25 = 8.2 * np.power(self.mass, 1 / 1.1) / 1e6
 
-        return PM_total
+        return (pm10, pm25)
diff --git a/dev/Input data_truck.xlsx b/dev/Input data_truck.xlsx
index bfe8a9b..38e9a2c 100644
Binary files a/dev/Input data_truck.xlsx and b/dev/Input data_truck.xlsx differ
diff --git a/setup.py b/setup.py
index 78e5864..420617e 100644
--- a/setup.py
+++ b/setup.py
@@ -28,7 +28,7 @@ def package_files(directory):
 
 setup(
     name="carculator_truck",
-    version="0.2.6",
+    version="0.2.7",
     packages=packages,
     author="Romain Sacchi <romain.sacchi@psi.ch>",
     license=open("LICENSE").read(),