{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Droplet distribution\n\nPlotting different droplet size distributions used in Opendrift\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom datetime import datetime, timedelta\nfrom opendrift.models.openoil import OpenOil"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "droplet size interval for plotting\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dmin = 1.e-6\ndmax = 3.e-3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "simulation with Johansen et al. (2015) droplet spectrum\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "o = OpenOil(loglevel=20, weathering_model='noaa')\no.set_config('environment:fallback:land_binary_mask', 0)\no.set_config('environment:fallback:x_sea_water_velocity', -.2)\no.set_config('environment:fallback:y_sea_water_velocity', 0)\no.set_config('environment:fallback:x_wind', 6.)\no.set_config('environment:fallback:y_wind', 0)\no.set_config('environment:fallback:sea_water_temperature', 5.)\no.set_config('environment:fallback:sea_surface_wave_stokes_drift_x_velocity', .3)\no.set_config('environment:fallback:sea_surface_wave_stokes_drift_y_velocity', 0)\no.set_config('wave_entrainment:droplet_size_distribution', 'Johansen et al. (2015)')\no.set_config('processes:evaporation', False)\no.set_config('processes:dispersion', False)\no.seed_elements(lon=4, lat=60, time=datetime.utcnow(), number=10000, radius=100,\n                 z=0, oil_type='TROLL, STATOIL', oil_film_thickness=0.005)\no.run(duration=timedelta(hours=2), time_step=3600)\n\ndroplet_diameters = o.elements.diameter\nsd = 0.4\nSd = np.log(10.)*sd\nDV_50 = np.median(droplet_diameters)\nDN_50 = np.exp( np.log(DV_50) - 3*Sd**2 )\n\nprint('DV_50: %f' % DV_50)\nprint('DN_50: %f' % DN_50)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=[14,14])\nplt.subplot(3, 2, 1)\nnVpdf, binsV, patches = plt.hist(droplet_diameters, 100, range=(dmin,dmax), align='mid')\nplt.xlabel('Droplet diameter d [m]', fontsize=8)\nplt.ylabel('V(d)', fontsize=8)\nplt.title('volume spectrum\\nJohansen et al. (2015) distribution in OpenOil', fontsize=10)\n\nplt.subplot(3,2,2)\nnVcum, binsV, patches = plt.hist(droplet_diameters, 100, range=(dmin,dmax), align='mid', cumulative=True)\nplt.xlabel('Droplet diameter d [m]', fontsize=8)\nplt.ylabel('V(d)', fontsize=8)\nplt.title('cumulative volume spectrum', fontsize=10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "calculate number spectrum from volume spectrum\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "d = 0.5* (binsV[1:] + binsV[:-1])\nV = 4./3. * np.pi * (d/2.)**3\nNpdf = nVpdf / V\nNcum = np.cumsum(Npdf)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "calculate theoretical cumulative Number distribution\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "spectrum = (np.exp(-(np.log(d) - np.log(DN_50))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from OpenOil median diameter\nDN_50_johansen = 0.000483\nspectrumJ = (np.exp(-(np.log(d) - np.log(DN_50_johansen))**2 / (2 * Sd**2))) / (d * Sd * np.sqrt(2 * np.pi)) # from parameters in Johansen et al. 2015\n\nspectrum_cum = np.cumsum(spectrum)\nspectrumJ_cum = np.cumsum(spectrumJ)\n# calculate theoretical number distribution pdf\n#spectrum_pdf = spectrum / np.sum(spectrum)\n#spectrumJ_pdf = spectrumJ / np.sum(spectrumJ)\n\n\nplt.subplot(3,2,3)\nplt.plot(d, Npdf/np.sum(Npdf), lw=2)\n#plt.plot(d, spectrum/np.sum(spectrum), label='curve fit from OpenOil', lw=2)\nplt.plot(d, spectrumJ/np.sum(spectrumJ), label='curve fit from Johansen et al. 2015', lw=2)\nplt.xlabel('Droplet diameter d [m]', fontsize=8)\nplt.ylabel('N(d)', fontsize=8)\nplt.title('number spectrum', fontsize=10)\n\nplt.subplot(3,2,4)\nplt.plot(d, Ncum/np.max(Ncum), label='OpenOil result', lw=2)\n#plt.plot(d, spectrum_cum/np.max(spectrum_cum), label='OpenOil par.', lw=2)\nplt.plot(d, spectrumJ_cum/np.max(spectrumJ_cum), label='Johansen et al. 2015', lw=2)\nplt.xlabel('Droplet diameter d [m]', fontsize=8)\nplt.ylabel('N(d)', fontsize=8)\nplt.title('cumulative number spectrum', fontsize=10)\nplt.legend(loc='lower right')\n\nplt.subplot(3,2,5)\nplt.loglog(d, Npdf/np.sum(Npdf), lw=2)\n#plt.loglog(d, spectrum/np.sum(spectrum), label='curve fit from OpenOil', lw=2)\nplt.loglog(d, spectrumJ/np.sum(spectrumJ), label='curve fit from Johansen et al. 2015', lw=2)\nplt.xlabel('Droplet diameter d [m]', fontsize=8)\nplt.ylabel('N(d)', fontsize=8)\nplt.title('number spectrum', fontsize=10)\n\nplt.subplot(3,2,6)\nplt.semilogx(d, Ncum/np.max(Ncum), label='OpenOil result', lw=2)\n#plt.semilogx(d, spectrum_cum/np.max(spectrum_cum), label='OpenOil par.', lw=2)\nplt.semilogx(d, spectrumJ_cum/np.max(spectrumJ_cum), label='Johansen et al. 2015', lw=2)\nplt.xlabel('Droplet diameter d [m]', fontsize=8)\nplt.ylabel('N(d)', fontsize=8)\nplt.title('cumulative number spectrum', fontsize=10)\n#plt.legend(loc='lower right')\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}