In [18]:
from atmPy import sizedistribution
from atmPy.instruments.DMA import smps
from atmPy.instruments.DMA import dma
import matplotlib.pyplot as plt
import atmPy.aerosol as aerosol
import atmPy.atmosphere as atm
import numpy as np
import pandas as pd
%matplotlib inline
The DMA is represented as a single object in the code. The object consists of three attributes: _ro, _ri and _l, the dimensions of the DMA (outer and inner radius - $r_o$ and $r_i$ as well as length $l$). These parameters are considered to be "private" and are defined by the children although the class can be used directly by setting these directly.
The DMA object contains one method: v2d. This method calculates the diameter of a size selected particle as a function of voltage, atmospheric conditions and flow rates.
def v2d(self, v, gas, qc, qm):
"""
Find selected diameter at a given voltage.
This function uses a Newton-Raphson root finder to solve for the
diameter of the particle.
Parameters
----------
v: float
Voltage in Volts
gas: gas object
Carrier gas object for performing calculations
qc: float
Input sheath flow in lpm
qm: float
Output sheath flow in lpm
Returns
-------
Diameter in nanometers
"""
gamma = self._l/log(self._ro/self._ri)
# Convert flow rates from lpm to m3/s
qc = float(qc)/60*0.001
qm = float(qm)/60*0.001
# Central mobility
zc = (qc+qm)/(4*pi*gamma*v)
return newton(lambda d: z(d, gas, 1)-zc, 1, maxiter=1000)
This method uses a Newton-Raphson zero finding routine to determine the diameter associated with central mobility.
The following code demonstrates how to use the DMA class that is found in the atmPy package.
In [21]:
# Instantiate a DMA object based on the NOAA wide DMA dimensions
noaa = dma.NoaaWide()
# We also need a gas object that will represent the carrier gas of interest: air
air = atm.Air()
# Set the conditions to something representative of the conditions in the Boulder DMA
air.t = 23
air.p = 820
d = noaa.v2d(50, air, 5, 5)
print('The diameter associated with a 5 lpm flow and 50 V is ' + str(d) + ' nm.')
In order to test the SMPS inversion, we need to simulate some type of size distribution
Generate a size distribution that we can use to test the retrieval of the SMPS data. We will put the geometric mean (dg) at 60 nm, close to the observed point during HAGiS.
In [2]:
sd = sizedistribution.simulate_sizedistribution(d=[10, 2000],
ndp=250,
dg=60,
sg=0.2,
nmode=3000)
f,a = sd.plot()
a.set_xlim([10,1000])
a.grid()
Now, we need to transform the size distribution so that it is what the SMPS might observe.
In [3]:
nd = sd.convert2numberconcentration()
f,a = nd.plot()
Now, we need to "correct" the size distribution so that it is what the SMPS would observe:
In [4]:
# demonstrate charing efficiency on size distribution
f2 = [aerosol.ndistr(i,n=-2,t=20) for i in nd.bincenters]
f3 = [aerosol.ndistr(i,n=-3,t=20) for i in nd.bincenters]
f1 = [aerosol.ndistr(i,n=-1,t=20) for i in nd.bincenters]
fig,ax = plt.subplots()
ax.plot(nd.bincenters, f1, 'r', nd.bincenters, f2,'k', nd.bincenters, f3, 'b')
ax.set_xscale('log')
ax.set_xlabel('Dp (nm)')
ax.set_ylabel('f')
ax.grid()
fig.tight_layout()
ax.set_title('Charging Efficiency vs. Diameter')
Out[4]:
In [5]:
# Create an air object of aerosl calculations...
air = atm.Air()
air.t = 20
air.p = 850
# Validate vs. Paul Barron's calculations; LOOKS GOOD
aerosol.z(1000, air, 1)
Out[5]:
In [6]:
# Define a function for finding the diameter
fmin = lambda dm: (np.abs(nd.bincenters - dm)).argmin()
xd = nd.copy()
for i,n in enumerate(xd.data.iloc[0,:].values):
f1 = aerosol.ndistr(xd.bincenters[i],n=-1,t = air.t)
f2 = aerosol.ndistr(xd.bincenters[i],n=-2,t = air.t)
f3 = aerosol.ndistr(xd.bincenters[i],n=-3,t = air.t)
d2 =aerosol.z2d(aerosol.z(xd.bincenters[i],air, n=2),air, n=1)
if d2>= xd.bins[0]:
k = fmin(d2)
if d2 > xd.bins[k]:
xd.data.values[0,k] += f2*n
else:
xd.data.values[0,k-1] +=f2*n
d3 =aerosol.z2d(aerosol.z(xd.bincenters[i],air, n=3),air, n=1)
if d3>= xd.bins[0]:
k = fmin(d3)
if d3 > xd.bins[k]:
xd.data.values[0,k] += f3*n
else:
xd.data.values[0,k-1] +=f3*n
xd.data.values[0,i] = f1*n
# print([i, xd.bincenters[i], fmin(d2), d2, fmin(d3), d3])
In [7]:
f,a = xd.plot()
a.grid()
So there is some funniness in the discreet nature of the problem - the reason the discreet peaks appear is that the algorithm is dumping multiples of multiples into a single bin. I have tried to address this by dumping the data into the bins with a fixed bin width (rather than searching for the nearest bin center), but to no avail. The problem is exacerbated by increasing the granularity of the problem (the opposite of what I thought might happen).
At this point, I will just move ahead and hope that I get this right...
In [8]:
# The SMPS routines need a mean data dataframe with the following entries
mean_data = {"Sh_Q_VLPM":3.9, "Aer_Temp_C": air.t, "Aer_Pres_PSI":air.p, "Aer_Q_VLPM":0.39
}
# Convert the dict to a list and generate a dataframe
df = pd.DataFrame([mean_data])
In [9]:
# Create a new SMPS object using the dimensions of the NOAA Wide DMA...
tSMPS = smps.SMPS(dma.NoaaWide())
In [10]:
xd.data.iloc[0,:].values
Out[10]:
In [11]:
dndlogdp = tSMPS.__fwhm__(xd.bincenters, xd.data.iloc[0,:].values, df)
In [12]:
dndlogdp[-1] = 0
In [13]:
xd.data.iloc[0,:].values
Out[13]:
In [14]:
xd.plot()
Out[14]:
In [15]:
fig, ax= plt.subplots()
ax.plot(tSMPS.diam_interp, dndlogdp)
ax.set_xscale('log')
ax.grid()
In [16]:
noaa = dma.NoaaWide()
air = atm.Air()
In [17]:
air.t = 25
air.p = 820
noaa.v2d(37.2, air, 3.7, 3.7)
Out[17]:
In [ ]: