In [1]:
import atmPy
from atmPy.aerosols.physics import hygroscopicity as hyg
from atmPy.data_archives import arm
from atmPy.aerosols.size_distribution import sizedistribution
%matplotlib inline
plt.style.use('hagen_default')

HygroscopicGrowthFactorDistributions

We need to generate a data set that can be used to initiate a HygroscopicGrowthFactorDistributions instance. Hiere we take Arm data generated by a HTDMA. The Arm data contains gf-distributions for different diameters, so we select one (200.0 nm).


In [2]:
fname = '../atmPy/unit_testing/test_data/sgptdmahygC1.b1.20120601.004227.cdf'
out = arm.read_netCDF(fname, data_quality= 'patchy', leave_cdf_open= False)

in general

Peaks in the gf-distribution are detected and fit by normal distributions (at log-base). Fit parameters are tightly constrained to avoid run-away parameters and unrealistic results, which in turn can result in unexpacted results ... hard coded fit parameters/boundaries might need adjustment.

Growth modes

Position of detected growthmodes and ratio of particles in it as a function of time. Here plotted on top of the gf-distribution time series.


In [3]:
out.hyg_distributions_d200nm.plot(growth_modes=True)


Out[3]:
(<matplotlib.figure.Figure at 0x110e14320>,
 <matplotlib.axes._subplots.AxesSubplot at 0x110e001d0>,
 <matplotlib.collections.QuadMesh at 0x10f70ffd0>,
 None)
/opt/local/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/colors.py:494: RuntimeWarning: invalid value encountered in less
  cbook._putmask(xa, xa < 0.0, -1)

In [5]:
out.hyg_distributions_d200nm.growth_modes_kappa


Out[5]:
ratio kappa
datetime
2012-05-31 23:57:27 0.032332 0.013393
2012-05-31 23:57:27 0.274847 0.093554
2012-05-31 23:57:27 0.692821 0.181959
2012-06-01 00:41:39 0.016535 0.004847
2012-06-01 00:41:39 0.349747 0.090270
2012-06-01 00:41:39 0.632345 0.187184
2012-06-01 00:41:39 0.001373 0.599974
2012-06-01 01:25:50 0.032276 0.017289
2012-06-01 01:25:50 0.160460 0.082738
2012-06-01 01:25:50 0.807264 0.187141
2012-06-01 02:10:05 0.039152 0.002397
2012-06-01 02:10:05 0.167292 0.074740
2012-06-01 02:10:05 0.726230 0.176295
2012-06-01 02:10:05 0.067326 0.224183
2012-06-01 02:54:20 0.060657 0.001958
2012-06-01 02:54:20 0.108566 0.058332
2012-06-01 02:54:20 0.830777 0.173118
2012-06-01 03:38:31 0.263576 0.005953
2012-06-01 03:38:31 0.542337 0.129795
2012-06-01 03:38:31 0.194087 0.183957
2012-06-01 04:22:45 0.133471 0.008650
2012-06-01 04:22:45 0.430195 0.104655
2012-06-01 04:22:45 0.436333 0.193954
2012-06-01 05:07:02 0.054511 0.004847
2012-06-01 05:07:02 0.011014 0.007807
2012-06-01 05:07:02 0.140996 0.074740
2012-06-01 05:07:02 0.793479 0.180480
2012-06-01 07:49:57 0.331540 0.029210
2012-06-01 07:49:57 0.668460 0.147077
2012-06-01 08:34:08 0.270126 0.017289
... ... ...
2012-06-01 15:12:02 0.588311 0.194143
2012-06-01 15:56:13 0.014403 -0.010766
2012-06-01 15:56:13 0.015722 -0.006278
2012-06-01 15:56:13 0.039883 0.029724
2012-06-01 15:56:13 0.929992 0.148641
2012-06-01 16:40:24 0.040341 0.002397
2012-06-01 16:40:24 0.959659 0.144723
2012-06-01 17:24:35 0.037704 0.005953
2012-06-01 17:24:35 0.789096 0.132626
2012-06-01 17:24:35 0.173199 0.225210
2012-06-01 18:08:46 0.052215 0.005121
2012-06-01 18:08:46 0.541978 0.116527
2012-06-01 18:08:46 0.392667 0.157896
2012-06-01 18:08:46 0.013140 0.415268
2012-06-01 18:52:57 0.068980 0.007807
2012-06-01 18:52:57 0.931020 0.127256
2012-06-01 19:37:07 0.115526 0.005953
2012-06-01 19:37:07 0.784259 0.129586
2012-06-01 19:37:07 0.100216 0.218123
2012-06-01 20:21:18 0.025811 -0.003602
2012-06-01 20:21:18 0.338658 0.102099
2012-06-01 20:21:18 0.635531 0.141822
2012-06-01 21:05:29 0.027640 0.005953
2012-06-01 21:05:29 0.246619 0.112243
2012-06-01 21:05:29 0.725740 0.157896
2012-06-01 21:49:40 0.037879 0.013393
2012-06-01 21:49:40 0.962121 0.137337
2012-06-01 22:33:52 0.070281 0.005484
2012-06-01 22:33:52 0.917147 0.111857
2012-06-01 22:33:52 0.012572 0.428707

86 rows × 2 columns

Mixing state

I came up with the following definition, it should be adjusted if there is a better one in the literature

Mixing state is given by the pythagoras of the particle ratios of all detected growth modes in a growth distribution. E.g. if there where three modes detected with ratios $r_1$, $r_2$, $r_3$ the mixing state is given by $\sqrt{r_1^2 + r_2^2 + r_3^2}$.


In [6]:
out.hyg_distributions_d200nm.mixing_state.plot(marker = 'o', ls = '')


Out[6]:
<matplotlib.axes._subplots.AxesSubplot at 0x110e39be0>

Grown size distribution

this is the sum. for optical properties the individual information is used so the change in the refractive index which is different for each growth mode is considered individually.


In [7]:
fname = '../atmPy/unit_testing/test_data/sgptdmaapssizeC1.c1.20120601.004227.cdf'
tdmaaps = arm.read_netCDF(fname, data_quality= 'patchy', leave_cdf_open= False)
sd = tdmaaps.size_distribution

In [8]:
hgfd = out.hyg_distributions_d200nm
# gmk = out.hyg_distributions_d200nm.growth_modes_kappa

In [9]:
sd.convert2dVdlogDp().plot()


Out[9]:
(<matplotlib.figure.Figure at 0x110bcee10>,
 <matplotlib.axes._subplots.AxesSubplot at 0x115172898>,
 <matplotlib.collections.QuadMesh at 0x115215d68>,
 <matplotlib.colorbar.Colorbar at 0x115434160>)

In [10]:
sd.hygroscopicity.parameters.growth_distribution = hgfd
sd.hygroscopicity.parameters.RH = 90

In [11]:
sd.hygroscopicity.grown_size_distribution.sum_of_all_sizeditributions.convert2dVdlogDp().plot()


Out[11]:
(<matplotlib.figure.Figure at 0x115e0dda0>,
 <matplotlib.axes._subplots.AxesSubplot at 0x115de6e80>,
 <matplotlib.collections.QuadMesh at 0x115df3dd8>,
 <matplotlib.colorbar.Colorbar at 0x115e57cf8>)

Optical properties

scattering


In [12]:
sd.optical_properties.parameters.wavelength = 550
sd.optical_properties.parameters.refractive_index = 1.5

In [13]:
sd.hygroscopicity.grown_size_distribution.optical_properties.scattering_coeff.plot()


/Users/htelg/prog/atm-py/atmPy/aerosols/physics/optical_properties.py:112: RuntimeWarning: invalid value encountered in true_divide
  y_phase_func = y_1p * 4 * _np.pi / scattering_cross_eff.sum()
Out[13]:
<matplotlib.axes._subplots.AxesSubplot at 0x116a10a90>

fRH


In [14]:
a = sd.hygroscopicity.f_RH_85_40.plot()
sd.hygroscopicity.f_RH_85_0.plot(ax = a)


/Users/htelg/prog/atm-py/atmPy/aerosols/physics/optical_properties.py:112: RuntimeWarning: invalid value encountered in true_divide
  y_phase_func = y_1p * 4 * _np.pi / scattering_cross_eff.sum()
Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x1170974a8>

catch fit runaways


In [11]:
from atmPy.tools import math_functions as _math_functions
def multi_gauss(x, *params, verbose=False):
    #     print(params)
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        if verbose:
            print(len(params), i)
        amp = params[i]
        pos = params[i + 1]
        sig = params[i + 2]
        y = y + _math_functions.gauss(x, amp, pos, sig)
    return y

In [ ]:
# %%debug --breakpoint /Users/htelg/prog/atm-py/atmPy/aerosols/physics/hygroscopicity.py:523
out.hyg_distributions_d200nm.plot(growth_modes=True)

atmPy.aerosols.physics.hygroscopicity ipdb> globals()['x'] = x ipdb> globals()['y'] = y ipdb> globals()['param'] = param ipdb> globals()['bound_l'] = bound_l ipdb> globals()['bound_h'] = bound_h globals()['x'] = x; globals()['y'] = y; globals()['param'] = param; globals()['bound_l'] = bound_l; globals()['bound_h'] = bound_h


In [7]:
x = atmPy.aerosols.physics.hygroscopicity.x
y = atmPy.aerosols.physics.hygroscopicity.y
param = atmPy.aerosols.physics.hygroscopicity.param
bound_l = atmPy.aerosols.physics.hygroscopicity.bound_l
bound_h= atmPy.aerosols.physics.hygroscopicity.bound_h

In [36]:
plt.plot(x,y)
plt.plot(x, y_start)
plt.plot(x, new_y)


Out[36]:
[<matplotlib.lines.Line2D at 0x1104aa780>]

In [ ]:
parry = np.array(param)
# parry[::3] *= 10

In [35]:
y_start = multi_gauss(x, *parry)

In [32]:
# fitres, _ = atmPy.aerosols.physics.hygroscopicity._curve_fit(multi_gauss, x, y, p0=param[:-3], bounds=(bound_l[:-3], bound_h[:-3]))
fitres, _ = atmPy.aerosols.physics.hygroscopicity._curve_fit(multi_gauss, x, y, p0=parry, bounds=(bound_l, bound_h), 
#                                                              max_nfev = 10000
                                                            )

new_y = multi_gauss(x, *fitres)

Kappa

In this section a kappa is defined instead of a growth distribution.


In [2]:
# fname = '../atmPy/unit_testing/test_data/sgptdmaapssizeC1.c1.20120601.004227.cdf
fname = '/Users/htelg/data/ARM/SGP/tdmaaps/sgptdmaapssizeC1.c1.20120201.002958.cdf'
tdmaaps = arm.read_netCDF(fname, data_quality= 'patchy', leave_cdf_open= False)

In [3]:
fname = '/Users/htelg/data/ARM/SGP/acsm/sgpaosacsmC1.b1.20120201.002022.cdf'
acsm = arm.read_netCDF(fname, data_quality= 'patchy', leave_cdf_open= False)

In [4]:
tdmaaps.size_distribution.parameters4reductions.wavelength = 550

In [5]:
# %%debug --breakpoint /Users/htelg/prog/atm-py/atmPy/aerosols/physics/hygroscopicity.py:606
tdmaaps.size_distribution.hygroscopicity.parameters.kappa = acsm.kappa
tdmaaps.size_distribution.hygroscopicity.parameters.refractive_index = acsm.refractive_index
fRH_nams_kams = tdmaaps.size_distribution.hygroscopicity.f_RH_85_0.copy()


/Users/htelg/prog/atm-py/atmPy/aerosols/instruments/AMS/AMS.py:63: UserWarning: Sulfate rich is not implemented yet. Mostly because I don't get it!
         There are not supposed to be any Nitrates or Chlorides present when we are in the sulfate rich regime ... but there are.
         I guess it has to do with the organics? Talk to chuck
  _warnings.warn(txt)
/Users/htelg/prog/atm-py/atmPy/aerosols/physics/hygroscopicity.py:52: RuntimeWarning: invalid value encountered in greater
  if _np.any(k > 1.4) or _np.any(k < 0):
/Users/htelg/prog/atm-py/atmPy/aerosols/physics/hygroscopicity.py:52: RuntimeWarning: invalid value encountered in less
  if _np.any(k > 1.4) or _np.any(k < 0):
/Users/htelg/prog/atm-py/atmPy/aerosols/physics/optical_properties.py:112: RuntimeWarning: invalid value encountered in true_divide
  y_phase_func = y_1p * 4 * _np.pi / scattering_cross_eff.sum()

In [6]:
# %%debug --breakpoint /Users/htelg/prog/atm-py/atmPy/aerosols/physics/hygroscopicity.py:742
tdmaaps.size_distribution.hygroscopicity.parameters.kappa = acsm.kappa
tdmaaps.size_distribution.hygroscopicity.parameters.refractive_index = 1.5#acsm.refractive_index
fRH_nfix_kams = tdmaaps.size_distribution.hygroscopicity.f_RH_85_0.copy()


/Users/htelg/prog/atm-py/atmPy/aerosols/physics/hygroscopicity.py:52: RuntimeWarning: invalid value encountered in greater
  if _np.any(k > 1.4) or _np.any(k < 0):
/Users/htelg/prog/atm-py/atmPy/aerosols/physics/hygroscopicity.py:52: RuntimeWarning: invalid value encountered in less
  if _np.any(k > 1.4) or _np.any(k < 0):
/Users/htelg/prog/atm-py/atmPy/aerosols/physics/optical_properties.py:112: RuntimeWarning: invalid value encountered in true_divide
  y_phase_func = y_1p * 4 * _np.pi / scattering_cross_eff.sum()

In [7]:
tdmaaps.size_distribution.hygroscopicity.parameters.kappa = acsm.kappa.data.values.mean() #acsm.kappa
tdmaaps.size_distribution.hygroscopicity.parameters.refractive_index = 1.5#acsm.refractive_index
fRH_nfix_kfix = tdmaaps.size_distribution.hygroscopicity.f_RH_85_0.copy()


/Users/htelg/prog/atm-py/atmPy/aerosols/physics/optical_properties.py:112: RuntimeWarning: invalid value encountered in true_divide
  y_phase_func = y_1p * 4 * _np.pi / scattering_cross_eff.sum()

In [8]:
a = fRH_nfix_kfix.plot(label = 'nfix kfix')
fRH_nfix_kams.plot(ax = a, label = 'nfix_kams')
fRH_nams_kams.plot(ax = a, label = 'nams_kams')
a.legend()


Out[8]:
<matplotlib.legend.Legend at 0x10b46a5c0>

In [ ]: