dust-bubble.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions\n",
    "\n",
    "|var name|value<br>exp1|value<br>exp2|unit|definition|\n",
    "|---|---|---|---|---|\n",
    "|posrlon|2.54517|same|[°E]||\n",
    "|posrlat|2.54879|same|[°N]|[coordinats](https://xkcd.com/2170/) of the emission point<br>(here in rotated coordinats of the forecast model)|\n",
    "|dt|30|same|[s]|model time step|\n",
    "|scale|2.e-7|same|-|scaling between measured particles and model particles<br>scale=2.e-7 1 model particle = 5 million measured particles|\n",
    "|emission time|10:50-11:45|13:40-14:30|-|time with particle emission|\n",
    "|ntmsp|110|100|-| number of time steps with particle emission|\n",
    "|itime|20.85|23.65|[h]| start time of the emission in hours after model start|\n",
    "|velox|0.|same|[ms-1]|the velocity of moving source in x direction|\n",
    "|veloy|0.|same|[ms-1]|the velocity of moving source in y direction|\n",
    "|expname|exp1|exp2|-|name of experiment|\n",
    "|plumeheight|5.|same|[m]|height of the dust plume|\n",
    "|plumeradius|5.|same|[m]|radius of the dust plume|\n",
    "|ngridpoints|100|same|-|number of horizontal gridpoints of the plume|\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posrlon=2.54517\n",
    "posrlat=2.54879\n",
    "\n",
    "dt=30\n",
    "\n",
    "scale=2.e-7\n",
    "\n",
    "velox=0.\n",
    "veloy=0.\n",
    "\n",
    "expname = 'exp1'\n",
    "ofilename=expname+'_'+str(scale)+'.out'\n",
    "\n",
    "# ntmsp\n",
    "if expname == 'exp1':\n",
    "    ntmsp = 110\n",
    "if expname == 'exp2':\n",
    "    ntmsp = 100\n",
    "\n",
    "# itime    \n",
    "if expname == 'exp1':\n",
    "    itime=20.85\n",
    "if expname == 'exp2':\n",
    "    itime=23.65\n",
    "    \n",
    "\n",
    "plumeheight=5.\n",
    "plumeradius=5.\n",
    "\n",
    "ngridpoints=100\n",
    "\n",
    "# filenames\n",
    "fname1 = 'GRIMM_1.5m_diff_'+expname+'.data'\n",
    "fname2 = 'GRIMM_3.8m_diff_'+expname+'.data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pre-calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def meter2degree(meter):\n",
    "    degree=np.zeros(len(meter))\n",
    "    for i in range(len(meter)):\n",
    "        degree[i] = meter[i]/6.3710088e6 * 180./np.pi\n",
    "    return degree\n",
    "\n",
    "\n",
    "dx=2.*plumeradius/ngridpoints # spacing of the gridboxes\n",
    "dz=dx \n",
    "\n",
    "# I divide the area form -plumeradius to +plumeradius in x and y direction \n",
    "# and from 0 to +plumeheight in z direction in grid boxes with the size dx**3\n",
    "ngridbox=int((ngridpoints+1)**2*(plumeheight/dz+1)) \n",
    "\n",
    "# list of the bubble midpoint position\n",
    "rlons=[posrlon+meter2degree([i*velox*dt])[0] for i in range(ntmsp)]\n",
    "rlats=[posrlat+meter2degree([i*veloy*dt])[0] for i in range(ntmsp)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input\n",
    "### Produce random data for testing\n",
    "The EDM data is not included in this repository but you can test the software with a randomized input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 data arrays with a len of 31\n",
    "\n",
    "dat1 = np.zeros(31)\n",
    "dat2 = np.zeros(31)\n",
    "\n",
    "npmax = 25000\n",
    "for i in range(len(dat1)):\n",
    "    dat1[i]=np.random.random_integers(0,npmax)    \n",
    "    dat2[i]=np.random.random_integers(0,npmax)\n",
    "    npmax -= npmax/5.\n",
    "\n",
    "binnames = ['0.25-0.28', '0.28-0.30', '0.30-0.35', '0.35-0.40', '0.40-0.45', '0.45-0.50', '0.50-0.58', \n",
    "            '0.58-0.65', '0.65-0.70', '0.70-0.80', '0.80-1.0', '1.0-1.3', '1.3-1.6', '1.6-2.0', '2.0-2.5', \n",
    "            '2.5-3.0', '3.0-3.5', '3.5-4.0', '4.0-5.0', '5.0-6.5', '6.5-7.5', '7.5-8.5', '8.5-10.0', \n",
    "            '10.0-12.5', '12.5-15.0', '15.0-17.5', '17.5-20.0', '20.0-25.0', '25.0-30.0', '30.0-32.0', '>32.0']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in data from Environmental Dust Monitor (EDM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# open input files\n",
    "f1 = open(fname1,encoding = 'unicode_escape')\n",
    "f2 = open(fname2,encoding = 'unicode_escape')\n",
    "\n",
    "# read header\n",
    "head=f1.readline()\n",
    "head=f1.readline()\n",
    "head=f2.readline()\n",
    "head=f2.readline()\n",
    "\n",
    "# read bin names form the header\n",
    "binnames = []\n",
    "for i in head.split():\n",
    "    if(int(i[0].isdigit())):\n",
    "        binnames.append(i)\n",
    "binnames.append('>32.0')\n",
    "\n",
    "# read the data to np.array\n",
    "data1 = np.loadtxt(f1,usecols=range(2,33))\n",
    "data2 = np.loadtxt(f2,usecols=range(2,33))\n",
    "\n",
    "#dimension of the data set\n",
    "dim_data=data1.shape\n",
    "\n",
    "\n",
    "# definition of the peak in the data & removal of the background noise\n",
    "\n",
    "# exp1\n",
    "if expname == 'exp1':   \n",
    "    peak1 = data1[454:464,:]-np.median(data1[454-50:464+50,:],axis=0)\n",
    "    peak2 = data2[458:465,:]-np.median(data2[458-50:465+50,:],axis=0)\n",
    "\n",
    "# exp2\n",
    "if expname == 'exp2':  \n",
    "    peak1 = data1[179:183,:]-np.median(data1[179-50:183+50,:],axis=0)\n",
    "    peak2 = data2[177:181,:]-np.median(data2[177-50:181+50,:],axis=0)\n",
    "\n",
    "\n",
    "peak1=np.where(peak1 > 0.,peak1,0.)\n",
    "peak2=np.where(peak2 > 0.,peak2,0.)\n",
    "\n",
    "dat1 = [np.sum(peak1[:,i])/(len(peak1)*6.) * dt for i in range(peak1.shape[1])]\n",
    "dat2 = [np.sum(peak2[:,i])/(len(peak2)*6.) * dt for i in range(peak2.shape[1])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the Bubble\n",
    "\n",
    "In this section, I use the EDM data to create a vertical profile of the particle concentration. \n",
    "From there I project the profile to a half-sphere to get a concentration field.\n",
    "Then I calculate how many particles there should be in every grid box of the sphere. \n",
    "Afterward, the particles are distributed randomly in the bubble. The probability that a particle appears at a certain position depends on the concentration field.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this lists store the data for the plotting later\n",
    "xp_plt=[]\n",
    "yp_plt=[]\n",
    "contour_plt=[]\n",
    "points_plt=[]\n",
    "nrmax_sv=[]\n",
    "ppgb_sv=[]\n",
    "nrsum=0 # sum of all particles in all bins \n",
    "\n",
    "\n",
    "# loop over all bins\n",
    "for ibin in range(len(dat1)):   # range(11,23) -> particles between 1 and 10 µm\n",
    "#    print(ibin)\n",
    "#    print(binnames[ibin])\n",
    "    \n",
    "    # negative particle concentrations may occur in some bins because of the background noise removal\n",
    "    # filter them out\n",
    "    if dat1[ibin] <= 0:\n",
    "        dat1[ibin] = 0.\n",
    "        \n",
    "    if dat2[ibin] <= 0:\n",
    "        dat2[ibin] = 0.\n",
    "        \n",
    "        \n",
    "    # define the data points we are using: surface:dat1, 1.5m:dat1, 3,8m:dat2, top of plume:0 \n",
    "    # it is defined as x points because we need a funktion y(x) for the fitting \n",
    "    # it is flipped later \n",
    "    xpoints=[dat1[ibin],dat1[ibin],dat2[ibin],0]\n",
    "    ypoints=[0.,1.5,3.8,plumeheight]\n",
    "\n",
    "\n",
    "\n",
    "    # fitting with 5th order Poly (5th order seems to work well)\n",
    "    x_p = np.linspace(0,plumeheight,int(plumeheight/dz)+1)\n",
    "    fit = np.polyfit(ypoints,xpoints,5)\n",
    "    y_p = np.polyval(fit,x_p)\n",
    "    \n",
    "    # filter out values <1 below the altitude with the peak value\n",
    "    maxpos=np.where(y_p == np.max(y_p))[0][0]\n",
    "    for i in range(maxpos):\n",
    "        if y_p[i] < 1.:\n",
    "            y_p[i] = 1.\n",
    "    \n",
    "\n",
    "    \n",
    "    # save data for plotting\n",
    "    xp_plt.append(x_p)\n",
    "    yp_plt.append(y_p)\n",
    "    \n",
    "    # the profile is now defined as funktion y_p(x_p). \n",
    "    # to make it a vertikal profile we flip it x_p(y_p)\n",
    "    \n",
    "    # not every profile reaches the possible plumeheight \n",
    "    # the next few lines search for the height where the profile becomes < 1\n",
    "    # this height is used later at the ellipsoid\n",
    "    zid=0\n",
    "    while y_p[zid] < 1.:\n",
    "        zid+=1\n",
    "    hmin=zid\n",
    "    while y_p[zid] >= 1.:\n",
    "        zid+=1\n",
    "        if zid == int(plumeheight/dz):\n",
    "            break\n",
    "    height=x_p[zid]\n",
    "\n",
    "    # the vertikal profile defines the particle conzentration \n",
    "    # in the center of the dust plume as a function of higth \n",
    "    \n",
    "    # now we use the profile to create a 3D conzentration field in the shape of a half sphere\n",
    "    # in the standart setting plumeradius=5. and plumehight=5.\n",
    "    # the conzentration field is defined in a grid from -5 to 5 meter horizontal and 0 to 5 meter vertical. \n",
    "    # with ngridpoints=100 each gridbox has a size of 0.1**3 meter\n",
    "    xg = np.linspace(-plumeradius, plumeradius,ngridpoints+1)\n",
    "    yg = np.linspace(-plumeradius, plumeradius,ngridpoints+1)\n",
    "    zg = np.linspace( 0.0, plumeheight,int(plumeheight/dz)+1)\n",
    "    \n",
    "\n",
    "    # get coordinates of the gridboxes with np.meshgrid\n",
    "    X, Y, Z = np.meshgrid(xg, yg,zg)\n",
    "\n",
    "    # I define the concentration field as an ellipsoid\n",
    "    # (x**2/a**2 + y**2/b**2 + z**2/c**2) = 1\n",
    "    # in this case, a and b is the plumeradius and c is the height\n",
    "    # This produces a normalized spherical concentration field with 1 in the centre and 0 at the edge\n",
    " \n",
    "    contour = 1 - (X**2/plumeradius**2  + Y**2/plumeradius**2  + Z**2/height**2 )\n",
    "    \n",
    "    # now we need to multiply the normalized concentration field with a vector \n",
    "    # that returns the vertical profile at the midpoint of the concentration field\n",
    "    # we search for the multiply value on every altitude level  \n",
    "    # multiply(level) = profiel value (level) / conture value (level) \n",
    "    \n",
    "    # multiply vector\n",
    "    multi=np.zeros(y_p.shape)\n",
    "    \n",
    "    # the midpoint of the bubble \n",
    "    mid=int(ngridpoints/2)\n",
    "    \n",
    "    # loop through level\n",
    "    for k in range(multi.shape[0]):\n",
    "        if contour[mid,mid,k] > 0.:\n",
    "            multi[k]=y_p[k]/contour[mid,mid,k]\n",
    "        else:\n",
    "            multi[k] = 0.\n",
    "            \n",
    "    # multiply the concentration field with the vector   \n",
    "    contour *= multi\n",
    "\n",
    "        \n",
    "    # if contour > 0. then contour == contour eles contour == 0.\n",
    "    contour = np.where(contour > 0., contour,0.)\n",
    "    \n",
    "    # remove possible artificials below the lower edge of the bubble\n",
    "    contour[:,:,:hmin]=0.0\n",
    "    \n",
    "    # save this contour field for plotting \n",
    "    contour_plt.append(contour)\n",
    "    \n",
    "    # ppgb - parts per grid box\n",
    "    # The conzentration [parts per liter] in transformed into parts per gridbox.\n",
    "    # size gridbox m**3 * part/liter * 1000 liter/m**3\n",
    "    \n",
    "    ppgb = dx**3 * contour*1000     \n",
    "    ppgb = ma.masked_where(ppgb < 0.,ppgb)\n",
    "    \n",
    "    # maximum number of particles in this bin. here the scaling took place \n",
    "    # (see table at the begin of this notebook)\n",
    "    # summ of particles * scale * number of emission time steps\n",
    "    nrmax=int(round(np.sum(ppgb)*scale)*ntmsp )   \n",
    "    \n",
    "    # save nrmax of this bin, this is needed later for the output \n",
    "    nrmax_sv.append(nrmax)\n",
    "    ppgb_sv.append(ppgb)\n",
    "    \n",
    "    # total sum of particles\n",
    "    nrsum+=nrmax\n",
    "    \n",
    "    # print particle number per bin\n",
    "    # print(nrmax)\n",
    "    \n",
    "    nrmax=nrmax_sv[ibin]\n",
    "    ppgb=ppgb_sv[ibin]\n",
    "    # 3d array that holds the coordinates of the paticles\n",
    "    points=np.zeros((nrmax,3))\n",
    "\n",
    "\n",
    "    # flatten the ppgb array and normalize it with the sum of particles\n",
    "    props=np.reshape(ppgb,ngridbox)/np.sum(ppgb)\n",
    "    npp=0\n",
    "    for pp in props.mask:\n",
    "        if pp:\n",
    "            props[npp]=0.0\n",
    "        npp+=1\n",
    "\n",
    "\n",
    "    # create nrmax random numbers in the range of ngridbox with the probabilty props\n",
    "    # So the single randmo number (ran in randoms) points to a specific cell in the flatten array\n",
    "    randoms = np.random.choice(ngridbox, nrmax,p=props)\n",
    "\n",
    "    nr=0\n",
    "    for ran in randoms:\n",
    "        # get the indices of the 3d array from the cell number of the flattened array \n",
    "        idx = np.unravel_index(np.ravel_multi_index((ran,), props.shape), ppgb.shape)\n",
    "        \n",
    "        # get the position of the point in cartesian coordinates \n",
    "        # then shift the particles randomly away from the middle point of the grid box\n",
    "        points[nr,0]=X[idx] + np.random.random()*dx-dx/2.\n",
    "        points[nr,1]=Y[idx] + np.random.random()*dx-dx/2.\n",
    "        points[nr,2]=Z[idx] + np.random.random()*dx-dx/2.\n",
    "        \n",
    "        # if the point are shifted underneath the surface, than get them back.\n",
    "        if points[nr,2] < 0.:\n",
    "            #print('warnig',points[nr,2])\n",
    "            points[nr,2]*=-1.\n",
    "        \n",
    "        \n",
    "        if points[nr,2] > 5.:\n",
    "            print('warnig',points[nr,2])\n",
    "            print(Z[idx],dx-dx/2.)\n",
    "\n",
    "            \n",
    "        nr+=1\n",
    "        \n",
    "    points_plt.append(points)\n",
    "    \n",
    "    \n",
    "print('total particles: ' + str(nrsum))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### particle numbers per bin\n",
    "\n",
    "|   bin   |exp1  |exp2  |\n",
    "|---------|------|------|\n",
    "|0.25-0.28|22440 |59500 | \n",
    "|0.28-0.30|24310 |34000 |\n",
    "|0.30-0.35|16390 |30000 |\n",
    "|0.35-0.40|7040  |29400 |\n",
    "|0.40-0.45|4730  |19500 |\n",
    "|0.45-0.50|3520  |7500  | \n",
    "|0.50-0.58|3410  |8800  | \n",
    "|0.58-0.65|2310  |7500  |\n",
    "|0.65-0.70|1100  |6600  |\n",
    "|0.70-0.80|3520  |12200 |\n",
    "|0.80-1.00|1320  |6600  |\n",
    "|1.00-1.30|1320  |6700  |\n",
    "|1.30-1.60|880   |3300  |\n",
    "|1.60-2.00|1100  |6500  |\n",
    "|2.00-2.50|1650  |13400 |\n",
    "|2.50-3.00|1430  |9700  |\n",
    "|3.00-3.50|440   |3000  |\n",
    "|3.50-4.00|440   |2100  |\n",
    "|4.00-5.00|880   |4200  |\n",
    "|5.00-6.50|660   |2500  |\n",
    "|6.50-7.50|440   |1100  |\n",
    "|7.50-8.50|220   |400   |\n",
    "|8.50-10.0|220   |700   |\n",
    "|10.0-12.5|220   |400   |\n",
    "|12.5-15.0|110   |300   |\n",
    "|15.0-17.5|110   |300   |\n",
    "|17.5-20.0|110   |200   |\n",
    "|20.0-25.0|110   |200   |\n",
    "|25.0-30.0|110   |200   |\n",
    "|30.0-32.0|110   |100   |\n",
    "|>32.0    |550   |300   |\n",
    "|total    |101200|277200|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Output \n",
    "\n",
    "The output file of this notebook is the start file of the Lagrangian particle dispersion model Itpas.\n",
    "The above-defined points are transformed into the rotated coordinates of the model.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def numspace(num):\n",
    "    nspace=0      # number of spaces between colums,\n",
    "                  # 3 for numbers < -10,\n",
    "                  # 4 for numbers < 0,\n",
    "                  # 5 for numbers >0,\n",
    "                  # 4 for numbers >10\n",
    "    if num <= -10.:\n",
    "        nspace=3\n",
    "        ndigs=8+4\n",
    "    elif num < 0.:\n",
    "        nspace=4\n",
    "        ndigs=8+3\n",
    "    elif num >= 0.:\n",
    "        nspace=5\n",
    "        ndigs=8+2\n",
    "        if num >= 10.:\n",
    "            nspace=4\n",
    "            ndigs=8+3\n",
    "    numlen=len(str(num))\n",
    "    dot=str(num).find('.')\n",
    "    afterdot=numlen-dot\n",
    "    nfill=0\n",
    "    if afterdot < 9:\n",
    "        nfill=9-afterdot\n",
    "    return nspace,ndigs,nfill\n",
    "\n",
    "\n",
    "\n",
    "# open new file\n",
    "f=open(ofilename,'w')\n",
    "\n",
    "# write header\n",
    "f.write(ofilename+' \\n')\n",
    "hline='time'\n",
    "hline+=5*' '+'lon'\n",
    "hline+=12*' '+'lat'\n",
    "hline+=12*' '+'height'\n",
    "hline+=7*' '+'npart'\n",
    "hline+=3*' '+'diam'\n",
    "hline+=14*' '+'dens'\n",
    "hline+=7*' '+'emission'\n",
    "hline+='\\n'\n",
    "f.write(hline)\n",
    "f.write(50*'-'+'\\n')\n",
    "\n",
    "\n",
    "# loop over all bins\n",
    "for ibin in range(len(dat1)):\n",
    "\n",
    "    # get points of the bin\n",
    "    points=points_plt[ibin]\n",
    "    \n",
    "    # get nrmax of the bin\n",
    "    nrmax=nrmax_sv[ibin]\n",
    "\n",
    "    # transform relative particle possitions to rotated coordinates\n",
    "    nrmax=int(nrmax/ntmsp)\n",
    "    for rlon,rlat,nr in zip(rlons,rlats,range(len(rlons))):\n",
    "        points[nr*nrmax:(nr+1)*nrmax,0]=rlon + meter2degree(points[nr*nrmax:(nr+1)*nrmax,0])\n",
    "        points[nr*nrmax:(nr+1)*nrmax,1]=rlat + meter2degree(points[nr*nrmax:(nr+1)*nrmax,1])\n",
    "\n",
    "\n",
    "    # define min and max particle diameter\n",
    "    pos=binnames[ibin].find('-')\n",
    "    if pos > 0:\n",
    "        Dmin=float(binnames[ibin][:pos])\n",
    "        Dmax=float(binnames[ibin][pos+1:])\n",
    "    else:\n",
    "        Dmin=float(binnames[ibin][1:])\n",
    "        Dmax=50.0\n",
    "        \n",
    "\n",
    "    # write points to ofile\n",
    "    npoint=0\n",
    "    time=0.0\n",
    "    initime=itime\n",
    "    for point in points:\n",
    "        npoint+=1\n",
    "        # create the output row\n",
    "        # time\n",
    "        if npoint> nrmax:\n",
    "            time+=0.0083\n",
    "            npoint=1\n",
    "\n",
    "        stime='{:8.4f}'.format(time+initime)\n",
    "        row=stime\n",
    "        # lon\n",
    "        nspace,ndigs,nfill=numspace(point[0])\n",
    "        row+=nspace*' '+str(point[0])[:ndigs]+nfill*'0'\n",
    "        #lat\n",
    "        nspace,ndigs,nfill=numspace(point[1])\n",
    "        row+=nspace*' '+str(point[1])[:ndigs]+nfill*'0'\n",
    "        # height\n",
    "        nspace,ndigs,nfill=numspace(point[2])\n",
    "        row+=nspace*' '+str(point[2])[:ndigs]+nfill*'0'\n",
    "        # npart\n",
    "        npart=1\n",
    "        nspace,ndigs,nfill=numspace(npart)\n",
    "        row+=nspace*' '+str(npart)\n",
    "        # Diam, chosen randomly in the size range of bin\n",
    "        diam=(np.random.random()*(Dmax-Dmin)+Dmin)*1e-6\n",
    "        nspace,ndigs,nfill=numspace(diam)\n",
    "        row+=nspace*' '+str(diam)[:ndigs]+str(diam)[-4:]\n",
    "        # Dens\n",
    "        dens=2650.0 #[kg/m**3]\n",
    "        nspace,ndigs,nfill=numspace(dens)\n",
    "        row+=nspace*' '+str(dens)\n",
    "        # emission\n",
    "        emission='0.00'\n",
    "        nspace,ndigs,nfill=numspace(float(emission))\n",
    "        row+=nspace*' '+emission\n",
    "        row+='\\n'\n",
    "        f.write(row)\n",
    "        \n",
    "# close file\n",
    "f.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "\n",
    "### Vertical profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open figure \n",
    "fig = plt.figure(figsize=(22,12)) \n",
    "\n",
    "# define grid of plots\n",
    "gs = gridspec.GridSpec(nrows=4, ncols=8)\n",
    "\n",
    "# x,y specify the subplot\n",
    "x=0\n",
    "y=0\n",
    "\n",
    "#loop over bins\n",
    "for ibin in range(len(dat1)):\n",
    "    \n",
    "    #open subplots\n",
    "    ax = fig.add_subplot(gs[y, x])\n",
    "\n",
    "    # set bin name as title\n",
    "    ax.set_title(binnames[ibin]+' $\\mu m$')\n",
    "    \n",
    "    # draw grid\n",
    "    ax.grid()\n",
    "\n",
    "    # get profile data\n",
    "    xpoints=[dat1[ibin],dat1[ibin],dat2[ibin],0]\n",
    "    ypoints=[0.,1.5,3.8,plumeheight]\n",
    "    \n",
    "    x_p=xp_plt[ibin]\n",
    "    y_p=yp_plt[ibin]\n",
    "    \n",
    "    # plot profile\n",
    "    ax.plot(xpoints,ypoints)\n",
    "    ax.plot(y_p,x_p)\n",
    "    ax.set_xlim(0)\n",
    "    ax.set_ylim(0)\n",
    "    \n",
    "\n",
    "    # increas x,y\n",
    "    x+=1\n",
    "    if x > 7:\n",
    "        x = 0\n",
    "        y+=1\n",
    "        \n",
    "fig.tight_layout()\n",
    "fig.savefig(ofilename[:-4]+'_profiles.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### Cross Section through the concentration field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open figure \n",
    "fig = plt.figure(figsize=(22,12))\n",
    "\n",
    "# define grid of plots\n",
    "gs = gridspec.GridSpec(nrows=4, ncols=8)\n",
    "\n",
    "# x,y specify the subplot\n",
    "x=0\n",
    "y=0\n",
    "\n",
    "# loop over bins\n",
    "for ibin in range(len(dat1)):\n",
    "    \n",
    "    #open subplots\n",
    "    ax = fig.add_subplot(gs[y, x])\n",
    "    \n",
    "    # set bin name as title\n",
    "    ax.set_title(binnames[ibin]+' $\\mu m$')\n",
    "    \n",
    "    # draw grid\n",
    "    ax.grid()\n",
    "\n",
    "    # get stored data\n",
    "    contour=contour_plt[ibin]\n",
    "    \n",
    "    # masked the areas where contour <= 0. so that they appear white in the plot\n",
    "    contour = ma.masked_where(contour <= 0.,contour)\n",
    "    \n",
    "    # the midpoint of the bubble \n",
    "    mid=int(ngridpoints/2-1)\n",
    "    # plot cross-section through the bubble\n",
    "    con = ax.contourf(X[0,:,:], Z[:,0,:], contour[:,mid,:],cmap='jet') #ppgb\n",
    "    fig.colorbar(con)\n",
    "\n",
    "\n",
    "    # increas x,y\n",
    "    x+=1\n",
    "    if x > 7:\n",
    "        x = 0\n",
    "        y+=1\n",
    "        \n",
    "fig.tight_layout()\n",
    "fig.savefig(ofilename[:-4]+'_contour.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3d scatter plot of the particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open figure \n",
    "fig = plt.figure(figsize=(22,12)) \n",
    "\n",
    "# define grid of plots\n",
    "gs = gridspec.GridSpec(nrows=4, ncols=8)\n",
    "\n",
    "# x,y specify the subplot\n",
    "x=0\n",
    "y=0\n",
    "\n",
    "#loop over bins\n",
    "for ibin in range(len(dat1)):\n",
    "    \n",
    "    #open subplots\n",
    "    ax = fig.add_subplot(gs[y, x],projection='3d')\n",
    "    \n",
    "    # set bin name as title\n",
    "    ax.set_title(binnames[ibin]+' $\\mu m$')\n",
    "    \n",
    "    # draw grid\n",
    "    ax.grid()\n",
    "    \n",
    "    # get points\n",
    "    points = points_plt[ibin]\n",
    "\n",
    "    # 3d scatter plot\n",
    "    ax.scatter(points[:,0],points[:,1],points[:,2],'b.',s=0.25)\n",
    "    #ax.set_xticks([])\n",
    "    #ax.set_yticks([])\n",
    "\n",
    "\n",
    "    # increas x,y\n",
    "    x+=1\n",
    "    if x > 7:\n",
    "        x = 0\n",
    "        y+=1\n",
    "        \n",
    "fig.tight_layout()\n",
    "fig.savefig(ofilename[:-4]+'_points.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### panel plot: Profile, cross-section, scatter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open figure \n",
    "fig = plt.figure(figsize=(12,4)) \n",
    "\n",
    "# define grid of plots\n",
    "gs = gridspec.GridSpec(nrows=1, ncols=3)\n",
    "\n",
    "# plot only bin 13 (1.6 - 2. µm)\n",
    "ibin = 10\n",
    "\n",
    "\n",
    "# open sub plots\n",
    "ax0 = fig.add_subplot(gs[0, 0])\n",
    "ax1 = fig.add_subplot(gs[0, 1])\n",
    "ax2 = fig.add_subplot(gs[0, 2],projection='3d')\n",
    "\n",
    "# ax0 profile\n",
    "x_p=xp_plt[ibin]\n",
    "y_p=yp_plt[ibin]\n",
    "\n",
    "xpoints=[dat1[ibin],dat1[ibin],dat2[ibin],0]\n",
    "ypoints=[0.,1.5,3.8,plumeheight]\n",
    "\n",
    "# plot profile\n",
    "\n",
    "ax0.plot(y_p,x_p)\n",
    "\n",
    "ax0.scatter(xpoints,ypoints,s=100)\n",
    "\n",
    "\n",
    "ax0.set_xlim(0)\n",
    "ax0.set_ylim(0,plumeheight+0.1)\n",
    "ax0.set_xlabel('particle per liter',size='x-large')\n",
    "ax0.set_ylabel('altitude [m]',size='x-large')\n",
    "\n",
    "ax0.annotate('Assumption #1',\n",
    "            xy=(xpoints[0], ypoints[0]), xycoords='data',\n",
    "            xytext=(-100, 25), textcoords='offset points',\n",
    "            arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3,rad=-.2\"))\n",
    "\n",
    "ax0.annotate('Assumption #2',\n",
    "            xy=(xpoints[3], ypoints[3]), xycoords='data',\n",
    "            xytext=(15, -35), textcoords='offset points',\n",
    "            arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3,rad=-.2\"))\n",
    "\n",
    "ax0.annotate('Measurements',\n",
    "            xy=(xpoints[1], ypoints[1]), xycoords='data',\n",
    "            xytext=(-80, 45), textcoords='offset points',\n",
    "            arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3,rad=.2\"))\n",
    "\n",
    "ax0.annotate('',\n",
    "            xy=(xpoints[2], ypoints[2]), xycoords='data',\n",
    "            xytext=(-73, -40), textcoords='offset points',\n",
    "            arrowprops=dict(arrowstyle=\"->\",connectionstyle=\"arc3,rad=-.3\"))\n",
    "\n",
    "# ax1 \n",
    "# get stored data\n",
    "contour=contour_plt[ibin]\n",
    "\n",
    "\n",
    "\n",
    "# masked the areas where contour <= 0. so that they appear white in the plot\n",
    "contour = ma.masked_where(contour <= 0.,contour)\n",
    "\n",
    "# the midpoint of the bubble \n",
    "mid=int(ngridpoints/2)\n",
    "# plot cross-section through the bubble\n",
    "con = ax1.contourf(X[0,:,:], Z[:,0,:], contour[:,mid,:],cmap='jet') #ppgb\n",
    "cb = fig.colorbar(con,ax=ax1)\n",
    "cb.ax.tick_params(labelsize='x-large')\n",
    "cb.set_label('particle per liter',size='x-large')\n",
    "\n",
    "ax1.set_ylim(0,plumeheight+0.1)\n",
    "ax1.set_ylabel('altitude [m]',size='x-large')\n",
    "ax1.set_xlabel('width [m]',size='x-large')\n",
    "\n",
    "# ax2\n",
    "# get points\n",
    "points = points_plt[ibin]\n",
    "\n",
    "# 3d scatter plot\n",
    "ax2.scatter(points[:,0],points[:,1],points[:,2],'b.',s=0.25)   \n",
    "#ax2.set_ylabel('width [m]',size='x-large')\n",
    "#ax2.set_zlabel('altitude [m]',size='x-large')\n",
    "#ax2.view_init(30, -65)\n",
    "\n",
    "ax0.tick_params(labelsize='x-large')\n",
    "ax1.tick_params(labelsize='x-large')\n",
    "ax2.tick_params(labelsize='x-large')\n",
    "\n",
    "ax0.set_title('a) Vertical profile',size='x-large')\n",
    "ax1.set_title('b) Concentration field',size='x-large')\n",
    "ax2.set_title('c) Particle positions',size='x-large')\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.savefig('out.png',dpi=500)"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}