From 146898fcf7b79269083644accbdb7ff35a445fe8 Mon Sep 17 00:00:00 2001
From: sfonxu <wronski.michal4@wp.pl>
Date: Fri, 18 Apr 2025 14:56:07 +0200
Subject: [PATCH] add first attempt at Twomey implementation to PySDM

---
 .../ccn_activation_spectrum/__init__.py       |  2 ++
 .../ccn_activation_spectrum/twomey_1959.py    | 15 ++++++++++
 .../physics/test_ccn_activation_spectrum.py   | 30 +++++++++++++++++++
 3 files changed, 47 insertions(+)
 create mode 100644 PySDM/physics/ccn_activation_spectrum/__init__.py
 create mode 100644 PySDM/physics/ccn_activation_spectrum/twomey_1959.py
 create mode 100644 tests/unit_tests/physics/test_ccn_activation_spectrum.py

diff --git a/PySDM/physics/ccn_activation_spectrum/__init__.py b/PySDM/physics/ccn_activation_spectrum/__init__.py
new file mode 100644
index 000000000..237fc2a9b
--- /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 000000000..5560fff7a
--- /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/tests/unit_tests/physics/test_ccn_activation_spectrum.py b/tests/unit_tests/physics/test_ccn_activation_spectrum.py
new file mode 100644
index 000000000..d68fb2367
--- /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