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')
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)
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.
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]:
In [5]:
out.hyg_distributions_d200nm.growth_modes_kappa
Out[5]:
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]:
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]:
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]:
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()
Out[13]:
In [14]:
a = sd.hygroscopicity.f_RH_85_40.plot()
sd.hygroscopicity.f_RH_85_0.plot(ax = a)
Out[14]:
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]:
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)
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()
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()
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()
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]:
In [ ]: