In [1]:
import sys,os,math
import numpy as np
import matplotlib as mpl
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cmx
#this works apparently only for savefig stuff
mpl.rcParams['figure.figsize']=(6.0,4.0) #(6.0,4.0)
mpl.rcParams['font.size']=10 #10
mpl.rcParams['savefig.dpi']=400 #72
mpl.rcParams['figure.subplot.bottom']=.1 #.125
plt.rc('font', family='serif')
plt.rc('text', usetex=True)
#inline Shit
%matplotlib inline
%config InlineBackend.figure_format='svg'
%config InlineBackend.rc = {'figure.facecolor': 'white', 'figure.subplot.bottom': 0.125, 'figure.edgecolor': 'white', 'savefig.dpi': 400, 'figure.figsize': (12.0, 8.0), 'font.size': 10}
#GUi shit
%matplotlib tk
mpl.get_configdir()
%load_ext autoreload
%autoreload 2
# Import General Packages from me
from Tools.Parsers import *
from Tools.BoundingBox import *
from Tools.Transformations import *
In [2]:
def fill_between_steps(x, y1, y2=0, h_align='mid', ax=None, **kwargs):
''' Fills a hole in matplotlib: Fill_between for step plots.
Parameters :
------------
x : array-like
Array/vector of index values. These are assumed to be equally-spaced.
If not, the result will probably look weird...
y1 : array-like
Array/vector of values to be filled under.
y2 : array-Like
Array/vector or bottom values for filled area. Default is 0.
**kwargs will be passed to the matplotlib fill_between() function.
'''
# If no Axes opject given, grab the current one:
if ax is None:
ax = plt.gca()
# First, duplicate the x values
xx = x.repeat(2)[1:]
# Now: the average x binwidth
xstep = (x[1:] - x[:-1]).mean()
# Now: add one step at end of row.
xx = np.append(xx, xx.max() + xstep)
# Make it possible to change step alignment.
if h_align == 'mid':
xx -= xstep / 2.
elif h_align == 'right':
xx -= xstep
# Also, duplicate each y coordinate in both arrays
y1 = y1.repeat(2)
if type(y2) == np.ndarray:
y2 = y2.repeat(2)
# now to the plotting part:
return ax.fill_between(xx, y1, y2=y2, **kwargs)
In [3]:
plt.close('all')
kR = np.array([1,5,20,80,97.5])/100;
sieveDiameter = np.array( [2.36,2,1.7,1.4,1.18]);
#linewitdh
lw = 2;
F = 1-kR;
print("kR:", kR)
print("F:", F)
print("sieveDiameter:", sieveDiameter)
sieveData = [ [x, (y - np.min(F)) / (np.max(F)-np.min(F)) ] for x,y in zip(sieveDiameter,F)]
sieveData = np.array(sorted(sieveData));
print("sieveData: \n", sieveData)
fm_d = np.diff(sieveData[:,1]) / np.diff(sieveData[:,0]);
print("density diameter: ", fm_d)
fig1, ax1 = plt.subplots(1,1, sharex=True)
ax1.grid(True)
ax1.set_ylim((0,3))
ax1.set_xlim((1,2.5))
ax1.set_title(r"diameter density")
ax1.set_xlabel("d [mm]")
ax1.set_ylabel(r" $\left[\frac{1}{\textnormal{mm}}\right]$")
x = np.hstack((sieveData[0,0],sieveData[:,0]));
y = np.hstack((0,fm_d,0));
ax1.plot(x,y, color='k', lw=lw, drawstyle='steps-post', label=r"$f_D^{m}$, $\mu = %.3f$, $d_{50}=1.554$" % mean_d_m)
fill_between_steps(sieveData[:,0],np.hstack((fm_d,0)),h_align='left', color=(0.4,0.4,0.4),alpha=0.5)
# Produce tranformed in number of particles
Fn_d_disk = np.zeros(np.shape(sieveData[:,0]));
Fn_d_disk[0] = 0;
# tranformiere mass diameter distribution zu anzahls disitribution
mean_d_n = 0;
mean_d_m = 0;
for i in range(0,len(f_d)):
print("generate transformed density for i=",i)
# integrate d * f_d_n over the whole range , f_d_n = d^-3*c_i \Xi_i(d) / normierung
mean_d_m += 0.5*fm_d[i]*( sieveData[i+1,0]**(2) - sieveData[i,0] **(2))
mean_d_n += fm_d[i]*(- sieveData[i+1,0]**(-1) + sieveData[i,0] **(-1))
# integrate f_d_n over the whole range , f_d_n = d^-3*c_i \Xi_i(d) / normierung
Fn_d_disk[i+1] = Fn_d_disk[i] + 0.5*fm_d[i] * ( - sieveData[i+1,0]**(-2) + sieveData[i,0] **(-2) )
mean_d_n /= Fn_d_disk[-1]
print("mean_d_m:", mean_d_m)
print("mean_d_n:", mean_d_n)
Fn_d_list = [];
fn_d_list = [];
for i in range(0,len(fm_d)):
d = np.linspace(sieveData[i,0],sieveData[i+1,0],100)
Fn_d_i = Fn_d_disk[i] + 0.5*fm_d[i] * ( - d**(-2) + sieveData[i,0] **(-2) )
fn_d_i = d**(-3) * fm_d[i]
# normierung mit F_d[-1]
Fn_d_list.append( (d, Fn_d_i / Fn_d_disk[-1]) );
fn_d_list.append( (d, fn_d_i / Fn_d_disk[-1]) );
Fn_d_disk = Fn_d_disk / Fn_d_disk[-1];
fig2, ax2 = plt.subplots(1,1, sharex=True)
ax2.margins(0.1)
ax2.grid(True)
ax2.set_title(r"cumulative distribution")
ax2.set_xlabel("d [mm]")
ax2.set_autoscaley_on(True)
ax2.plot(sieveData[:,0],sieveData[:,1],'k-', label=r"$F_D^{m}$" ,lw=lw )
ax2.plot(sieveData[:,0],sieveData[:,1],'ko' )
d_total = np.array([sieveData[0,0]]);
fn_total = np.array([0]);
for d,f in fn_d_list:
fn_total = np.hstack((fn_total,f))
d_total = np.hstack((d_total,d))
d_total = np.hstack((d_total,[sieveData[-1,0]]));
fn_total = np.hstack((fn_total,[0]));
ax1.plot(d_total,fn_total,'b',lw=lw, label=r"$f_D^{n}$, $\mu = %.3f$, $d_{50}=1.481$" % mean_d_n)
ax1.fill_between(d_total, 0, fn_total, facecolor='blue', alpha=0.2)
d_total = np.array(());
Fn_total = np.array(());
for d,F in Fn_d_list:
d_total = np.hstack((d_total,d))
Fn_total = np.hstack((Fn_total,F))
ax2.plot(d_total,Fn_total,'b-', label=r'$F_D^{n}$', lw=lw)
ax2.plot(sieveData[:,0],Fn_d_disk,'bo')
ax2.plot(sieveData[:,0],Fn_d_disk,'b--', label=r'$\bar{F}_D^{n}$')
nrRange = [(i,sieveData[i,0],f) for (i,),f in np.ndenumerate(Fn_d_disk) ]
for x in nrRange[1:-1]:
ax2.annotate('(%s,%1.3f)' % (x[1],x[2]) , xy=x[1:3], textcoords='offset points', horizontalalignment='right', verticalalignment='bottom') # <--
fn_d_disk = np.diff(Fn_d_disk) / np.diff(sieveData[:,0]);
median_d_n_disk = (0.5-Fn_d_disk[1])*(sieveData[2,0]-sieveData[1,0])/(Fn_d_disk[2]-Fn_d_disk[1]) + sieveData[1,0];
rho_for_sim = 1/(median_d_n_disk*1e-3**3*math.pi/6)
print("rho_for_sim:",rho_for_sim)
mean_d_n_disk=0;
for i in range(0,len(fn_d_disk)):
mean_d_n_disk += 0.5*fn_d_disk[i]*( sieveData[i+1,0]**(2) - sieveData[i,0] **(2))
print("fn_d_disk:",fn_d_disk)
x = np.hstack((sieveData[0,0],sieveData[:,0]));
y = np.hstack((0,fn_d_disk,0));
ax1.plot(x,y,'b--', drawstyle='steps-post', label=r"$\bar{f}_D^{m}$, $\mu=%.3f$, $d_{50}=%.3f$" % (mean_d_n_disk,median_d_n_disk))
leg1 = ax1.legend(loc='upper right')
leg2 = ax2.legend(loc='lower right')
for legobj in leg1.legendHandles + leg2.legendHandles:
legobj.set_linewidth(lw)
fig1.tight_layout()
fig2.tight_layout()
fig1.canvas.draw()
fig2.canvas.draw()
fig1.savefig("densStarlitebead.pdf")
fig2.savefig("cumdistStarlitebead.pdf")
In [ ]:
2**(-2)
In [ ]:
sieveData[:-1,0]
In [ ]:
Fn_d_disk[:-1]
In [ ]:
Fn_d_disk
In [ ]: