In [131]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Dp = 4.
rhow = 1025.
rhos = 2650.
nf = 2.
nd = 15
D = np.logspace(np.log10(20.),np.log10(1500.),nd)
# density from fractal dimension
rhof = rhow + (rhos-rhow)*(D/Dp)**(nf-3.)
# log-normal distribution
mu = np.log10(600.)
sigma = np.log10(2.)
fr = 1./(np.sqrt(2.*np.pi*sigma**2))*np.exp(-(np.log10(D)-mu)**2/(2.*sigma**2))
fr = fr/(np.sum(fr))
# de-floc distribution
frd = np.zeros_like(fr)
frd[0:5]=(.1, .2, .4, .2, .1)
for t in np.arange(0.,1.1,.1):
frt = fr+t*(frd-fr)
print t, np.sum(frt), frt
for i in np.arange(len(D)):
print D[i],rhof[i],fr[i]
fig = plt.figure(figsize=(8,10))
plt.subplot(2,1,1)
plt.plot(D,fr)
plt.subplot(2,1,2)
plt.plot(D,rhof)
plt.show()
In [129]:
nDbins = 21
nrbins = 31
Dbinlims = np.logspace(np.log10(2.),np.log10(2000.),nDbins+1)
rbinlims = np.linspace(1030.,1350.,nrbins+1)
# print rbinlims
# print Dbinlims
frac = np.zeros( (nDbins,nrbins) )
for k in np.arange(len(D)):
for i in np.arange(len(Dbinlims)-1):
for j in np.arange(len(rbinlims)-1):
if (D[k]>Dbinlims[i] and D[k]<=Dbinlims[i+1]) and (rhof[k]>rbinlims[j] and rhof[k]<=rbinlims[j+1]):
print k,i,j,D[k],rhof[k],frt[k],'rbinlims:',rbinlims[j],rbinlims[j+1]
frac[i,j]=frac[i,j]+fr[k]
print 'Sum of frac: ',np.sum(frac)
frac = frac/np.sum(frac)
In [120]:
Dbins = Dbinlims[0]+np.cumsum(0.5*np.diff(Dbinlims))
rbins = rbinlims[0]+np.cumsum(0.5*np.diff(rbinlims))
print Dbins
print rbins
x,y = np.meshgrid( rbins, Dbins )
print x.shape, y.shape, frac.shape
x = x.flatten()
y = y.flatten()
z = frac.flatten()
import matplotlib.cm as cm
import matplotlib.colors as colors
norm = colors.Normalize(frac.min(), frac.max())
colors = cm.hsv(norm(z))
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection='3d')
ax.bar3d( x, y, np.zeros(len(z)), 1, 20, z, color = colors)
plt.ylabel('Size (microns)')
plt.xlabel('Density (kg/m^3)')
#plt.zlabel('Fraction')
plt.show()