diff --git a/PySDM/formulae.py b/PySDM/formulae.py
index 93040d286a..57afd142c9 100644
--- a/PySDM/formulae.py
+++ b/PySDM/formulae.py
@@ -57,6 +57,7 @@ def __init__(  # pylint: disable=too-many-locals
         terminal_velocity: str = "GunnKinzer1949",
         air_dynamic_viscosity: str = "ZografosEtAl1987",
         bulk_phase_partitioning: str = "Null",
+        ccn_activation_spectrum: str = "Null",
         handle_all_breakups: bool = False,
     ):
         # initialisation of the fields below is just to silence pylint and to enable code hints
@@ -90,7 +91,8 @@ def __init__(  # pylint: disable=too-many-locals
         self.air_dynamic_viscosity = air_dynamic_viscosity
         self.terminal_velocity = terminal_velocity
         self.bulk_phase_partitioning = bulk_phase_partitioning
-
+        self.ccn_activation_spectrum = ccn_activation_spectrum
+        
         self._components = tuple(
             i
             for i in dir(self)
diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py
index 65ad336cbe..b429cbde72 100644
--- a/PySDM/physics/__init__.py
+++ b/PySDM/physics/__init__.py
@@ -47,5 +47,6 @@
     air_dynamic_viscosity,
     terminal_velocity,
     bulk_phase_partitioning,
+    ccn_activation_spectrum,
 )
 from .constants import convert_to, in_unit, si
diff --git a/PySDM/physics/ccn_activation_spectrum/__init__.py b/PySDM/physics/ccn_activation_spectrum/__init__.py
new file mode 100644
index 0000000000..237fc2a9b0
--- /dev/null
+++ b/PySDM/physics/ccn_activation_spectrum/__init__.py
@@ -0,0 +1,2 @@
+from PySDM.impl.null_physics_class import Null
+from .twomey_1959 import Twomey1959
diff --git a/PySDM/physics/ccn_activation_spectrum/twomey_1959.py b/PySDM/physics/ccn_activation_spectrum/twomey_1959.py
new file mode 100644
index 0000000000..5560fff7ac
--- /dev/null
+++ b/PySDM/physics/ccn_activation_spectrum/twomey_1959.py
@@ -0,0 +1,15 @@
+import numpy as np
+"""
+[Twomey 1959](https://doi.org/10.1007/BF01993560). Note that supersaturation is expressed in temperature unit as the elevation of the dew point or as percent of supersaturation; and concentrations are reported for 10C and 800 mb. 
+"""
+
+class Twomey1959:
+    def __init__(self, const):
+        assert np.isfinite(const.TWOMEY_K)
+        assert np.isfinite(const.TWOMEY_N0)
+
+    @staticmethod
+    def ccn_concentration(const, saturation_ratio):
+        return const.TWOMEY_N0*np.power(saturation_ratio-1, const.TWOMEY_K)
+    
+    
diff --git a/PySDM/physics/constants_defaults.py b/PySDM/physics/constants_defaults.py
index 5f4a19a85b..2802e9f5f9 100644
--- a/PySDM/physics/constants_defaults.py
+++ b/PySDM/physics/constants_defaults.py
@@ -684,3 +684,34 @@ def compute_derived_values(c: dict):
     c["water_molar_volume"] = c["Mv"] / c["rho_w"]
     c["rho_STP"] = c["p_STP"] / c["Rd"] / c["T_STP"]
     c["H_u"] = c["M"] / c["p_STP"]
+
+W76W_G0 = -2.9912729e3 * si.K**2
+W76W_G1 = -6.0170128e3 * si.K
+W76W_G2 = 1.887643854e1
+W76W_G3 = -2.8354721e-2 * si.K**-1
+W76W_G4 = 1.7838301e-5 * si.K**-2
+W76W_G5 = -8.4150417e-10 * si.K**-3
+W76W_G6 = 4.4412543e-13 * si.K**-4
+W76W_G7 = 2.858487
+W76W_G8 = 1 * si.Pa
+
+B80W_G0 = 6.112 * si.hPa
+B80W_G1 = 17.67 * si.dimensionless
+B80W_G2 = 243.5 * si.K
+
+one_kelvin = 1 * si.K
+
+bulk_phase_partitioning_T_cold = 235 * si.K
+bulk_phase_partitioning_T_warm = 273 * si.K
+bulk_phase_partitioning_exponent = np.nan
+
+
+BOLIN_ISOTOPE_TIMESCALE_COEFF_C1 = np.nan * si.dimensionless
+"""
+Coeffitient c1 used in [Bolin 1958](https://https://digitallibrary.un.org/record/3892725)
+for the falling drop evaporation timescale of equilibration with ambient air void of a given
+isotopologue; in the paper timescale is calculated for tritium with assumption of no tritium
+in the environment around the drop (Table 1).
+"""
+TWOMEY_K = np.nan
+TWOMEY_N0 = np.nan
diff --git a/tests/unit_tests/physics/test_ccn_activation_spectrum.py b/tests/unit_tests/physics/test_ccn_activation_spectrum.py
new file mode 100644
index 0000000000..d68fb23671
--- /dev/null
+++ b/tests/unit_tests/physics/test_ccn_activation_spectrum.py
@@ -0,0 +1,30 @@
+import numpy as np
+from matplotlib import pyplot
+
+from PySDM import Formulae
+from PySDM.physics.constants import PER_CENT, si, in_unit
+
+def test_twomey_and_wojciechowski_1969_fig1(plot=True):
+    """[Twomey, Wojciechowski 1969](https://doi.org/10.1175/1520-0469(1969)26%3C648:OOTGVO%3E2.0.CO;2)
+    """
+    #arange
+    for k, N in zip([0.5,0.7], [100,500]):
+        formulae = Formulae(ccn_activation_spectrum="Twomey1959",constants={"TWOMEY_K": k, "TWOMEY_N0": N/si.cm**3})
+        supersaturation = np.logspace(np.log10(.2), np.log10(9))*PER_CENT
+        #act
+        activated_nuclei_concentration = formulae.ccn_activation_spectrum.ccn_concentration(saturation_ratio=supersaturation+1)
+        #plot
+        pyplot.plot(in_unit(supersaturation,PER_CENT),
+                    in_unit(activated_nuclei_concentration,si.cm**-3),
+                    label=f"{k=}"
+                )
+    pyplot.xlim(0.1,10)
+    pyplot.ylim(1,1000)
+    pyplot.xscale("log")
+    pyplot.yscale("log")
+    pyplot.xlabel("Percent supersaturation")
+    pyplot.ylabel("Nuclei [cm$^{-3}$]")
+    pyplot.grid()
+    pyplot.legend()
+    pyplot.show()
+    #assert