In [3]:
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))
# 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 [19]:
nDbins = 15
nrbins = 7
Dbinlims = np.logspace(np.log10(4.),np.log10(1800.),nDbins+1)
rbinlims = np.linspace(1030.,1350.,nrbins+1)
print rbinlims
rbinlims = np.array((1025.,1030.,1050.,1080.,1240.,1300.,1500.,2700.))
print rbinlims
# print Dbinlims
# de-floc distribution
frd = np.zeros_like(fr)
frd[0:5]=(.1, .2, .4, .2, .1)
t = .0
fn = "defloc_t{0}.png".format(100.*t)
print fn
frt = fr+t*(frd-fr)
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]+frt[k]
print 'Sum of frac: ',np.sum(frac)
frac = frac/np.sum(frac)
# plot histogram
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.afmhot(norm(z))
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection='3d')
ax.bar3d( x, y, np.zeros(len(z)), 20, 30, z, color = colors)
plt.ylabel('Size (microns)',fontsize=14)
plt.xlabel('Density (kg/m^3)',fontsize=14)
plt.title('t={0}'.format(t))
ax.set_zlim([0.,.45])
#plt.zlabel('Fraction')
plt.savefig(fn)
plt.show()
In [6]:
# test relaxation calcs
frd = np.zeros_like(fr)
frd[0:5]=(.1, .2, .4, .2, .1)
for t in np.arange(0., 1.2, .2):
fr = fr+t*(frd-fr)
print fr
In [ ]: