From 9195e67a3ba6e541fc4257b2a2ecfcae27616449 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Mon, 14 Jul 2025 21:34:58 +0200 Subject: [PATCH 1/2] fractionation factors for different S and temp --- .../Zaba_et_al/fractionation_factors.ipynb | 169 ++++++++++++++++++ 1 file changed, 169 insertions(+) create mode 100644 examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb diff --git a/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb b/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb new file mode 100644 index 000000000..fefc5150a --- /dev/null +++ b/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb @@ -0,0 +1,169 @@ +{ + "cells": [ + { + "cell_type": "code", + "id": "initial_id", + "metadata": { + "collapsed": true, + "ExecuteTime": { + "end_time": "2025-07-14T18:42:23.295115Z", + "start_time": "2025-07-14T18:42:23.287266Z" + } + }, + "source": [ + "import numpy as np\n", + "from matplotlib import pyplot\n", + "from functools import partial\n", + "from open_atmos_jupyter_utils import show_plot\n", + "\n", + "from PySDM import Formulae\n", + "from PySDM.physics import si" + ], + "outputs": [], + "execution_count": 71 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-14T19:16:41.463007Z", + "start_time": "2025-07-14T19:16:41.449038Z" + } + }, + "cell_type": "code", + "source": [ + "formulae = Formulae(\n", + " isotope_equilibrium_fractionation_factors=\"MerlivatAndNief1967+Majoube1970+Majoube1971\",\n", + " isotope_kinetic_fractionation_factors=\"JouzelAndMerlivat1984\",\n", + " isotope_diffusivity_ratios=\"Stewart1975\"\n", + ")\n", + "\n", + "saturation = np.linspace(.7, 1.1, 7)\n", + "temperature_C = np.linspace(-30, 0)\n", + "temperature = formulae.trivia.C2K(temperature_C) * si.K\n" + ], + "id": "5f8f76c80c1d2286", + "outputs": [], + "execution_count": 126 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-14T19:16:41.836941Z", + "start_time": "2025-07-14T19:16:41.833451Z" + } + }, + "cell_type": "code", + "source": [ + "alpha_liquid = formulae.isotope_equilibrium_fractionation_factors.alpha_l_2H(temperature)\n", + "alpha_kinetic = partial(formulae.isotope_kinetic_fractionation_factors.alpha_kinetic, alpha_equilibrium=alpha, D_ratio_heavy_to_light=formulae.isotope_diffusivity_ratios.ratio_2H_heavy_to_light(temperature))\n", + "alpha_ice = formulae.isotope_equilibrium_fractionation_factors.alpha_i_2H(temperature)" + ], + "id": "810fcbb7b502850d", + "outputs": [], + "execution_count": 127 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-14T19:34:08.839335Z", + "start_time": "2025-07-14T19:34:07.644009Z" + } + }, + "cell_type": "code", + "source": [ + "for S in saturation:\n", + " alpha = alpha_liquid * alpha_kinetic(saturation=S)\n", + " pyplot.plot(alpha, temperature_C, 'k', alpha=0.1)\n", + " pyplot.annotate(f'{S:.3g}', xy=(alpha[4], temperature_C[4]), rotation=55, color='k', alpha=0.5, size=10)\n", + "\n", + "pyplot.plot(alpha_kinetic(saturation=1)*alpha_ice, temperature_C, label='wrt ice')\n", + "print(alpha[3]/alpha[0])\n", + "pyplot.plot(alpha_kinetic(saturation=1)*alpha_liquid, temperature_C, label='wrt liquid')\n", + "pyplot.gca().set(\n", + " xlabel='Fractionation factor',\n", + " ylabel=\"Temperature [K]\",\n", + ")\n", + "pyplot.gca().invert_yaxis()\n", + "pyplot.legend()\n", + "show_plot()\n" + ], + "id": "2f1aeb9076b01e6f", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.9965711674381411\n" + ] + }, + { + "data": { + "text/plain": [ + "
" + ], + "image/svg+xml": "\n\n\n \n \n \n \n 2025-07-14T21:34:08.818040\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "HBox(children=(HTML(value=\"./tmp_l8_pbgl.pdf
\"), HTML(value…" + ], + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "a8f58c2a3af844b08fb503e47f0af6f2" + } + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "execution_count": 149 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-14T19:23:49.618983Z", + "start_time": "2025-07-14T19:23:49.615602Z" + } + }, + "cell_type": "code", + "source": "", + "id": "3a1b25d03fa7c846", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": "", + "id": "79642bd18579086a" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 93a965f425b5836a8ddb8557f2ddef86ffec26a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Agnieszka=20=C5=BBaba?= Date: Tue, 15 Jul 2025 21:06:08 +0200 Subject: [PATCH 2/2] fix args naming --- .../Zaba_et_al/fractionation_factors.ipynb | 76 +++++++++---------- 1 file changed, 36 insertions(+), 40 deletions(-) diff --git a/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb b/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb index fefc5150a..7193ed0cb 100644 --- a/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb +++ b/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb @@ -6,8 +6,8 @@ "metadata": { "collapsed": true, "ExecuteTime": { - "end_time": "2025-07-14T18:42:23.295115Z", - "start_time": "2025-07-14T18:42:23.287266Z" + "end_time": "2025-07-15T18:53:32.653518Z", + "start_time": "2025-07-15T18:53:30.861075Z" } }, "source": [ @@ -20,13 +20,13 @@ "from PySDM.physics import si" ], "outputs": [], - "execution_count": 71 + "execution_count": 1 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-07-14T19:16:41.463007Z", - "start_time": "2025-07-14T19:16:41.449038Z" + "end_time": "2025-07-15T18:53:33.014638Z", + "start_time": "2025-07-15T18:53:32.737631Z" } }, "cell_type": "code", @@ -38,47 +38,58 @@ ")\n", "\n", "saturation = np.linspace(.7, 1.1, 7)\n", + "\n", + "# temperature_C might be replaced with temperature profile\n", "temperature_C = np.linspace(-30, 0)\n", "temperature = formulae.trivia.C2K(temperature_C) * si.K\n" ], "id": "5f8f76c80c1d2286", "outputs": [], - "execution_count": 126 + "execution_count": 2 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-07-14T19:16:41.836941Z", - "start_time": "2025-07-14T19:16:41.833451Z" + "end_time": "2025-07-15T18:53:34.101574Z", + "start_time": "2025-07-15T18:53:33.022297Z" } }, "cell_type": "code", "source": [ - "alpha_liquid = formulae.isotope_equilibrium_fractionation_factors.alpha_l_2H(temperature)\n", - "alpha_kinetic = partial(formulae.isotope_kinetic_fractionation_factors.alpha_kinetic, alpha_equilibrium=alpha, D_ratio_heavy_to_light=formulae.isotope_diffusivity_ratios.ratio_2H_heavy_to_light(temperature))\n", - "alpha_ice = formulae.isotope_equilibrium_fractionation_factors.alpha_i_2H(temperature)" + "alpha_l = formulae.isotope_equilibrium_fractionation_factors.alpha_l_2H(temperature)\n", + "alpha_l_kinetic = partial(\n", + " formulae.isotope_kinetic_fractionation_factors.alpha_kinetic,\n", + " alpha_equilibrium=alpha_l,\n", + " D_ratio_heavy_to_light=formulae.isotope_diffusivity_ratios.ratio_2H_heavy_to_light(temperature)\n", + ")\n", + "\n", + "alpha_i = formulae.isotope_equilibrium_fractionation_factors.alpha_i_2H(temperature)\n", + "alpha_i_kinetic = partial(\n", + " formulae.isotope_kinetic_fractionation_factors.alpha_kinetic,\n", + " alpha_equilibrium=alpha_i,\n", + " D_ratio_heavy_to_light=formulae.isotope_diffusivity_ratios.ratio_2H_heavy_to_light(temperature)\n", + ")" ], "id": "810fcbb7b502850d", "outputs": [], - "execution_count": 127 + "execution_count": 3 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-07-14T19:34:08.839335Z", - "start_time": "2025-07-14T19:34:07.644009Z" + "end_time": "2025-07-15T18:53:34.418783Z", + "start_time": "2025-07-15T18:53:34.110259Z" } }, "cell_type": "code", "source": [ "for S in saturation:\n", - " alpha = alpha_liquid * alpha_kinetic(saturation=S)\n", + " alpha = alpha_l * alpha_l_kinetic(saturation=S)\n", " pyplot.plot(alpha, temperature_C, 'k', alpha=0.1)\n", " pyplot.annotate(f'{S:.3g}', xy=(alpha[4], temperature_C[4]), rotation=55, color='k', alpha=0.5, size=10)\n", "\n", - "pyplot.plot(alpha_kinetic(saturation=1)*alpha_ice, temperature_C, label='wrt ice')\n", - "print(alpha[3]/alpha[0])\n", - "pyplot.plot(alpha_kinetic(saturation=1)*alpha_liquid, temperature_C, label='wrt liquid')\n", + "pyplot.plot(alpha_i_kinetic(saturation=1)*alpha_i, temperature_C, label='wrt ice')\n", + "pyplot.plot(alpha_l_kinetic(saturation=1)*alpha_l, temperature_C, label='wrt liquid')\n", "pyplot.gca().set(\n", " xlabel='Fractionation factor',\n", " ylabel=\"Temperature [K]\",\n", @@ -89,19 +100,12 @@ ], "id": "2f1aeb9076b01e6f", "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "0.9965711674381411\n" - ] - }, { "data": { "text/plain": [ "
" ], - "image/svg+xml": "\n\n\n \n \n \n \n 2025-07-14T21:34:08.818040\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" + "image/svg+xml": "\n\n\n \n \n \n \n 2025-07-15T20:53:34.400056\n image/svg+xml\n \n \n Matplotlib v3.10.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n" }, "metadata": {}, "output_type": "display_data" @@ -109,40 +113,32 @@ { "data": { "text/plain": [ - "HBox(children=(HTML(value=\"./tmp_l8_pbgl.pdf
\"), HTML(value…" + "HBox(children=(HTML(value=\"./tmphss6ofj6.pdf
\"), HTML(value…" ], "application/vnd.jupyter.widget-view+json": { "version_major": 2, "version_minor": 0, - "model_id": "a8f58c2a3af844b08fb503e47f0af6f2" + "model_id": "99cfa7e0f64b411abb28c15e1bea5957" } }, "metadata": {}, "output_type": "display_data" } ], - "execution_count": 149 + "execution_count": 4 }, { "metadata": { "ExecuteTime": { - "end_time": "2025-07-14T19:23:49.618983Z", - "start_time": "2025-07-14T19:23:49.615602Z" + "end_time": "2025-07-15T18:53:34.461114Z", + "start_time": "2025-07-15T18:53:34.459575Z" } }, "cell_type": "code", "source": "", - "id": "3a1b25d03fa7c846", + "id": "79642bd18579086a", "outputs": [], "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "outputs": [], - "execution_count": null, - "source": "", - "id": "79642bd18579086a" } ], "metadata": {