Trying to mimic a mathematical model based on the following model documentation.

4.1.1.4 Exposure to Vapour: Evaporation

https://www.rivm.nl/bibliotheek/rapporten/2017-0197.pdf

Python script: 

import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

class ExposureToVapourEvaporation:
    def __init__(self):
        # Input parameters
        self.frequency = 197  # per year
        self.exposure_duration = 0.75  # minute
        self.product_amount = 500  # g (amount of diluted product applied on a surface)
        self.weight_fraction_substance = 0.2  # fraction in the product
        self.room_volume = 1  # m³
        self.ventilation_rate = 0.5  # per hour
        self.inhalation_rate = 22.9 / 1000  # m³/min (converted from L/min)
        self.vapour_pressure = 0.0106  # Pa
        self.molecular_weight = 46.1 / 1000  # kg/mol (converted from g/mol)
        self.release_area = 0.002  # m²
        self.release_duration = 0.3  # minute
        self.application_temperature = 20 + 273.15  # K (converted from °C)
        self.mass_transfer_coefficient = 10 / 3600  # m/s (converted from m/h)
        self.body_weight = 9.8  # kg
        
        # Optional parameters
        self.is_product_used_in_dilution = True
        self.dilution = 1  # times (only used if is_product_used_in_dilution is True)
        self.molecular_weight_matrix = 22 / 1000  # kg/mol (converted from g/mol)
        self.is_pure_substance = False
        
        # Derived parameters
        self.is_constant_surface_area = True
        self.weight_fraction_solution = self.calculate_weight_fraction_solution()

    def calculate_weight_fraction_solution(self):
        if self.is_product_used_in_dilution:
            return self.weight_fraction_substance / self.dilution
        return self.weight_fraction_substance

    def calculate_equilibrium_vapour_pressure(self):
        return self.vapour_pressure

    def evaporation_ode(self, y, t):
        A_air, A_prod = y
        K = self.mass_transfer_coefficient
        P_eq = self.calculate_equilibrium_vapour_pressure()
        P_air = A_air * 8.314 * self.application_temperature / (self.molecular_weight * self.room_volume)
        
        if t < self.release_duration * 60:  # During release
            dA_air_dt = K * self.release_area * (self.molecular_weight / (8.314 * self.application_temperature)) * (P_eq - P_air) - (self.ventilation_rate / 3600) * self.room_volume * A_air
            dA_prod_dt = -dA_air_dt

        elif t == self.release_duration * 60:  # At the exact end of release
            dA_air_dt = 100 * K * self.release_area * (self.molecular_weight / (8.314 * self.application_temperature)) * (P_eq - P_air) - (self.ventilation_rate / 3600) * self.room_volume * A_air
        

        
        else:  # After release
            dA_air_dt = -(self.ventilation_rate / 3600) * self.room_volume * A_air
            dA_prod_dt = 0
        
        if A_prod <= 0:
            dA_prod_dt = 0
            dA_air_dt = -(self.ventilation_rate / 3600) * self.room_volume * A_air
        
        return [dA_air_dt, dA_prod_dt]


    def solve_evaporation(self):
        t = np.linspace(0, self.exposure_duration * 60, 1000)  # seconds
        t = np.sort(np.unique(np.append(t, self.release_duration * 60)))  # Ensure we have a point at the end of release
        y0 = [0, (self.product_amount / 1000) * self.weight_fraction_solution]  # kg
        solution = odeint(self.evaporation_ode, y0, t)
        return t, solution

    def calculate_metrics(self):
        t, solution = self.solve_evaporation()
        A_air = solution[:, 0]
        
        concentrations = A_air / self.room_volume * 1e6  # Convert to mg/m³
        mean_concentration = np.mean(concentrations)
        peak_concentration = np.max(concentrations)
        
        # Calculate TWA 15 min
        twa_15_min = mean_concentration if self.exposure_duration <= 15 else np.mean(concentrations[-int(15*60/self.exposure_duration):])
        
        # Calculate daily and yearly averages
        daily_average = mean_concentration * (self.exposure_duration / 1440)
        yearly_average = daily_average * (self.frequency / 365)
        
        # Calculate doses
        event_dose = (mean_concentration * self.inhalation_rate * self.exposure_duration) / self.body_weight
        day_dose = event_dose  # Assuming one event per day
        
        metrics = {
            "mean_event_concentration": mean_concentration,
            "peak_concentration_twa_15_min": twa_15_min,
            "mean_concentration_day": daily_average,
            "year_average_concentration": yearly_average,
            "external_event_dose": event_dose,
            "external_day_dose": day_dose
        }
        
        return metrics, t, concentrations

    def plot_air_concentration(self, t, concentrations):
        plt.figure(figsize=(10, 6))
        plt.plot(t / 60, concentrations)  # Convert time to minutes
        plt.axvline(x=self.release_duration, color='r', linestyle='--', label='End of Release')
        plt.axvline(x=self.exposure_duration, color='g', linestyle='--', label='End of Exposure')
        plt.title('Air Concentration Over Time')
        plt.xlabel('Time (minutes)')
        plt.ylabel('Concentration (mg/m³)')
        plt.legend()
        plt.grid(True)
        plt.yscale('log')
        plt.show()


# Create an instance of the model and calculate metrics
model = ExposureToVapourEvaporation()
metrics, t, concentrations = model.calculate_metrics()

print("Calculated Results:")
for key, value in metrics.items():
    print(f"{key}: {value:.2e} mg/m³" if "concentration" in key else f"{key}: {value:.2e} mg/kg bw")

# Plot the air concentration over time
model.plot_air_concentration(t, concentrations)

The expected result is:

Current:

So the Spike should occur at the point of the End of release.

Been trying for ages to get the maths to work but just can’t get the jump to match at 0.3.