myprojects-hg-reborn

Fork

#! /usr/bin/env python
from scipy import *
import pylab  # used to read the .csv file

slip_yannis=1


#T_0=273.+150.

save_knu=1

x=pylab.load('V_list') # x stands for the volume list
delta_V=x[2]-x[1]

D_f=pylab.load('fractal_dimension') # value of the fractal dimension
print "the fractal dimension is", D_f

a_0=pylab.load('monomer_radius')


diam_seq=(pylab.load('D_mean_seq'))/1e9 # I also load the diameter grid
# which I convert into m since I am using SI units



lensum=len(x) # number of particle bins
print 'lensum is', lensum

#Now a list of physical parameters or constants which do not change according to the expe
#rimental conditions
rho_p=1.*pylab.load('rho_0') # density [kg/m^3] of the particles	
R=0.05 # anaconda radius in m
k_B=1.38e-23 # Boltzmann const in SI units



air_properties=pylab.load("output_air") # I load some air properties calculated
#  at fixed T with another code. Useful above all when I work at T=const inside the
# anaconda


T_0=air_properties[0]
print "the air temperature is [K]", T_0
lam_air=air_properties[1] # value of air mean-free path [in m]. 
print "lam_air is, ", lam_air

#Treated as a fixed parameter right now
nu_air=air_properties[2] #kinematic air viscosity
print 'the air kinematic viscosity is [m/s]', nu_air
mu=air_properties[3] # air dynamic viscosity
print 'the air dynamic viscosity is [Kg/(m**2*s)]', mu



def Brow_ker_fuchs_class_optim_slip(Vlist,diam_seq):
	knu=2.*lam_air/diam_seq # vector with the Knudsen numbers
	#print 'knu is', knu
        if (save_knu==1):
            pylab.save("knudsen",knu)
	if (slip_yannis==0):
		Diff=k_B*T_0/(3.*pi*mu*diam_seq)
	elif (slip_yannis==1):
		#print 'k_B is ', k_B
		#print 'T is, ', T_0
		Diff=k_B*T_0/(3.*pi*mu*diam_seq)*(1.+knu*(1.17+0.53*exp(-0.78/knu)))
	#print 'Diff is', Diff
	m=rho_p*Vlist # this holds in general (i.e. for Vlist !=3.)
	## as long as Vlist is the volume of solid
	##and not the space occupied by the agglomerate structure
	c=(8.*k_B*T_0/(pi*m))**0.5
	#print 'c is', c
	l=8.*Diff/(pi*c)
        if (save_knu==1):
            pylab.save("lambda_part",l)
	#print 'l is', l
	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq

	beta=(\
	(diam_seq[:,newaxis]+diam_seq[newaxis,:])   \
	/(diam_seq[:,newaxis]+diam_seq[newaxis,:]\
	+2.*(g[:,newaxis]**2.+g[newaxis,:]**2.)**0.5)\
	+ 8.*(Diff[:,newaxis]+Diff[newaxis,:])/\
	((c[:,newaxis]**2.+c[newaxis,:]**2.)**0.5*\
	(diam_seq[:,newaxis]+diam_seq[newaxis,:]))\
	)**(-1.)
	
	#print 'the mean of beta is', mean(beta)
	
	## now I have all the bits for the kernel matrix
	kern_mat=Brow_ker_cont_optim_slip(Vlist,diam_seq)*beta

	#print 'beta is', beta


	
	return kern_mat


def Brow_ker_cont_optim_slip(Vlist,diam_seq):
	# same expression for the kernel in the continuous limit
	# as the one used by Maricq.
	knu=2.*lam_air/diam_seq # vector with the Knudsen numbers
	if (slip_yannis ==0):
		Slip=zeros(len(knu))
		Slip[:]=1.
	elif (slip_yannis==1):
		Slip=1.+knu*(1.17+0.53*exp(-0.78/knu))
	kern_mat=2.*k_B*T_0/(3.*mu)*(Vlist[:,newaxis]**(1./D_f)+\
	Vlist[newaxis,:]**(1./D_f))* \
	(Slip[:,newaxis]/Vlist[:,newaxis]**(1./D_f)+\
	Slip[newaxis,:]/Vlist[newaxis,:]**(1./D_f))
	return kern_mat

def Brow_ker_free_optim(Vlist):
	lam=2./D_f-0.5
	kern_mat=(3./(4.*pi))**lam*(6.*k_B*T_0/rho_p)**0.5*a_0**(2.-6./D_f)*\
	(1./Vlist[:,newaxis]+1./Vlist[newaxis,:])**0.5*\
	(Vlist[:,newaxis]**(1./D_f)+Vlist[newaxis,:]**(1./D_f))**2.
	return kern_mat




def beta_calc(Vlist,diam_seq):
	knu=2.*lam_air/diam_seq # vector with the Knudsen numbers
	#print 'knu is', knu
	if (slip_yannis==0):
		Diff=k_B*T_0/(3.*pi*mu*diam_seq)*((5.+4.*knu+6.*knu**2.+18.*knu**3.)/\
		(5.-knu+(8.+pi)*knu**2.))
	elif (slip_yannis==1):
		#print 'k_B is ', k_B
		#print 'T is, ', T_0
		Diff=k_B*T_0/(3.*pi*mu*diam_seq)*(1.+knu*(1.17+0.53*exp(-0.78/knu)))
	#print 'Diff is', Diff
	m=rho_p*Vlist # this holds in general (i.e. for Vlist !=3.)
	## as long as Vlist is the volume of solid
	##and not the space occupied by the agglomerate structure
	c=(8.*k_B*T_0/(pi*m))**0.5
	#print 'c is', c
	l=8.*Diff/(pi*c)
	#print 'l is', l
	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq

	beta=(\
	(diam_seq[:,newaxis]+diam_seq[newaxis,:])   \
	/(diam_seq[:,newaxis]+diam_seq[newaxis,:]\
	+2.*(g[:,newaxis]**2.+g[newaxis,:]**2.)**0.5)\
	+ 8.*(Diff[:,newaxis]+Diff[newaxis,:])/\
	((c[:,newaxis]**2.+c[newaxis,:]**2.)**0.5*\
	(diam_seq[:,newaxis]+diam_seq[newaxis,:]))\
	)**(-1.)
	
	#print 'the mean of beta is', mean(beta)
	
	## now I have all the bits for the kernel matrix
	

	#print 'beta is', beta
	

	
	return beta




print "1st kernel calculation"



my_kernel=Brow_ker_fuchs_class_optim_slip(x,diam_seq)

print "after 1st kernel calculation"


#It is convenient to carry out the integration in log space

#Now I want to create the matrix M_mat[i,j]=K[i,j]n[i]n[j]

n=pylab.load('population_binned') # Now I load the initial state of the system
N_tot=n.sum()
print "the total population is, ",N_tot

M_mat=my_kernel*n[:,newaxis]*n[newaxis,:]
M_mat2=M_mat

#Now it is time to evaluate the matrix N_mat[i,j]=n[i]n[j]
N_mat=n[:,newaxis]*n[newaxis,:]
N_mat2=N_mat
#mmmmhhhhh....I want to carry out the integration in log space; since my binning in
#diameters is lognormal

M_mat=M_mat*diam_seq[:,newaxis]*diam_seq[newaxis,:]
N_mat=N_mat*diam_seq[:,newaxis]*diam_seq[newaxis,:]


#Now I carry out the mean:

Mean_ker_on_dist=M_mat.sum(axis=0).sum(axis=0)/N_mat.sum(axis=0).sum(axis=0)


Mean_ker_on_dist2=M_mat2.sum(axis=0).sum(axis=0)/N_mat2.sum(axis=0).sum(axis=0)


print "the mean kernel averaged on the initial distr is, ",Mean_ker_on_dist
print "the mean kernel averaged on the initial distr (correct!) is, ",Mean_ker_on_dist2


mean_part=pylab.load("mean_particle")
print 'mean_part is', mean_part

#print "the kernel evaluated on the mean particle is"

mean_dia_vec=zeros(len(diam_seq))
mean_vol_vec=zeros(len(diam_seq))


mean_dia_vec[:]=mean_part[1]
mean_vol_vec[:]=mean_part[0]


save_knu=0
			
my_aver_ker=Brow_ker_fuchs_class_optim_slip(mean_vol_vec,mean_dia_vec)

#print "the kernel evaluated for the average particle is, ", my_aver_ker[0,0]

tau_mean=1./(Mean_ker_on_dist*sum(n))
tau_mean2=1./(Mean_ker_on_dist2*sum(n))


print "the distribution-averaged tau is, ", tau_mean


print "the CORRECT distribution-averaged tau is, ", tau_mean2


#Now a calculation of the beta factor

beta_fac= beta_calc(x,diam_seq)

#print "the shape of beta is, ", shape(beta_fac)

#print "beta is, ", beta_fac
pylab.loglog(diam_seq*1e9,diag(beta_fac),linewidth=2.)
#xlim=(diam_seq.min()*1e10,diam_seq.max()*1e9)
pylab.xlabel('diameter [nm]')
pylab.ylabel('beta factor')
#pylab.legend(('prey population','predator population'))
pylab.title('beta vs diameter')
pylab.grid(True)
pylab.savefig("beta_evolution.pdf")
pylab.hold(False)
pylab.clf()

#Now I calculate the free molecular kernel


free_kernel=Brow_ker_free_optim(x)

cont_kernel=Brow_ker_cont_optim_slip(x,diam_seq)

slip_yannis=0

cont_kernel_no_slip=Brow_ker_cont_optim_slip(x,diam_seq)

print "the free-kernel for the smallest particle is, ", free_kernel[0,0]
print "and Fuchs for the smallest particle is, ", my_kernel[0,0]

my_knu=pylab.load("knudsen")

# mark_knu=where(abs(my_knu-1.)==min(abs(my_knu-1.)),diam_seq,0.)

# select_diam=max(mark_knu)


# print "select_diam is,", select_diam*1e9


pylab.loglog(diam_seq*1e9,my_knu)
pylab.xlabel('diameter [nm]')
pylab.ylabel('knudsen number')
#pylab.legend(('prey population','predator population'))
pylab.title('knudsen number vs diameter')
pylab.grid(True)
pylab.savefig('test_knudsen.pdf')
pylab.hold(False)
pylab.clf()


#Now I plot of the diagonal along the kernel

pylab.loglog(diam_seq*1e9,diag(my_kernel),diam_seq*1e9,diag(free_kernel)\
             ,diam_seq*1e9,diag(cont_kernel),diam_seq*1e9,diag(cont_kernel_no_slip))
#pylab.axvline(x=select_diam*1e9)
pylab.xlabel('diameter [nm]')
pylab.ylabel('kernel value along diagonal')
pylab.legend(('Fuchs','Free Molecular','Continuum','continuum no slip'))
pylab.title('Kernel number vs diameter')
pylab.grid(True)
pylab.savefig('test_kernel_diagonal.pdf')
pylab.hold(False)
pylab.clf()

my_l=pylab.load("lambda_part")

lam_air_vec=zeros(len(my_l))
lam_air_vec[:]=lam_air*1e9
		  

pylab.loglog(diam_seq*1e9,my_l*1e9,diam_seq*1e9,lam_air_vec)
pylab.xlabel('diameter [nm]')
pylab.ylabel('lambda particle')
#pylab.legend(('prey population','predator population'))
pylab.title('Lambda Brownian particle vs diameter')
pylab.grid(True)
pylab.savefig('test_lambda_part.pdf')
pylab.hold(False)
pylab.clf()


print "So far so good"