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..7193ed0cb --- /dev/null +++ b/examples/PySDM_examples/Zaba_et_al/fractionation_factors.ipynb @@ -0,0 +1,165 @@ +{ + "cells": [ + { + "cell_type": "code", + "id": "initial_id", + "metadata": { + "collapsed": true, + "ExecuteTime": { + "end_time": "2025-07-15T18:53:32.653518Z", + "start_time": "2025-07-15T18:53:30.861075Z" + } + }, + "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": 1 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-15T18:53:33.014638Z", + "start_time": "2025-07-15T18:53:32.737631Z" + } + }, + "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", + "\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": 2 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-15T18:53:34.101574Z", + "start_time": "2025-07-15T18:53:33.022297Z" + } + }, + "cell_type": "code", + "source": [ + "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": 3 + }, + { + "metadata": { + "ExecuteTime": { + "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_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_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", + ")\n", + "pyplot.gca().invert_yaxis()\n", + "pyplot.legend()\n", + "show_plot()\n" + ], + "id": "2f1aeb9076b01e6f", + "outputs": [ + { + "data": { + "text/plain": [ + "
" + ], + "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" + }, + { + "data": { + "text/plain": [ + "HBox(children=(HTML(value=\"./tmphss6ofj6.pdf
\"), HTML(value…" + ], + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "99cfa7e0f64b411abb28c15e1bea5957" + } + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "execution_count": 4 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-15T18:53:34.461114Z", + "start_time": "2025-07-15T18:53:34.459575Z" + } + }, + "cell_type": "code", + "source": "", + "id": "79642bd18579086a", + "outputs": [], + "execution_count": null + } + ], + "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 +}