myprojects-hg-reborn

Fork

#! /usr/bin/env python
import scipy as s
import pylab as p  # used to read the .csv file
import scipy.integrate as si
from scipy import stats #I need this module for the linear fit

#x=pylab.load('V_list') # x stands for the volume list


r_mon=0.5

v_mono=4./3.*s.pi*r_mon**3.

print "the monodisperse volume is, ", v_mono


n_mon=s.logspace(s.log10(1.), s.log10(2500.),200)   # volume grid

x=n_mon*v_mono #volume of solids expressed in terms of number of monomers




print "n_mon is, ", n_mon

m=n_mon #since each monomer has mass 1 in my units

threshold=15. # I define the threshold between small and large clusters, to be used where it matters.


small=s.where(n_mon<=threshold)
large=s.where(n_mon>threshold)



a_small=0.36
df_small=2.20

a_large=0.24
df_large=1.62


a_tot=0.278

df_tot=1.72



        






n_part=5000.  #total number of monomers in the box

rho=0.01  #monomer density in the box

choice_kernel= 10  #0: standard continuum kernel (and standard diffusion_coefficient)
                  #1: constant kernel
                  #2: continuum kernel, 2 df's and special treatment of monomers (special diffusion coeff)
                  #3: continuum kernel, 1 df and special treatment of monomers (special diffusion coeff)
                  #4: continuum kernel, 2 df's without special treatment of
                  #monomers (special diffusion coeff)
                  #5: continuum kernel, 2 df's, Fuchs beta without special treatment of monomers
                  #(special diffusion coeff)
                  #6: continuum kernel, 1 df, Fuchs beta without special treatment of monomers
                  #(special diffusion coeff)
                  #7: continuum kernel, 2 df's, Fuchs beta without special treatment of monomers
                  #(standard diffusion coeff)
                  #8: continuum kernel, 1 df, Fuchs beta without special treatment of monomers
                  #(standard diffusion coeff)
                  #9: continuum kernel, 2 df's, Fuchs beta without special treatment of monomers, but
                  # Fuchs correction is applied to monomers only.
                  #10: as in 5, but this time we introduce by hand the radius of gyration of the smallest
                  #clusters
                  


choice_slope = 2 #2: 2 different df's
                 #1: 1 different df


t=s.linspace(0.,3000.,3000) # I am choosing my time grid for time evolution




t_fin=len(t)-1



beta=1. #cluster/monomer 1/tau
k_B=1. #in these units
T_0=0.5 #temperature of the system
m_mon=1. #monomer mass in these units
sigma=1. #monomer diameter
mu=(m_mon*beta)/(3.*s.pi*sigma) # fluid viscosity





# End of the setting up of the code. Now I use the objects I created above


########################################################################################
########################################################################################
########################################################################################
########################################################################################
########################################################################################
########################################################################################



def fitting_stat(x,y):
    slope, intercept, r, prob2, see = stats.linregress(x,y)
    #print "len(x) is, ", len(x)
    if (len(x)>2):
        
        see=see*s.sqrt(len(x)/(len(x)-2.))

        mx = x.mean()
        sx2 = ((x-mx)**2).sum()
        sd_intercept = see * s.sqrt(1./len(x) + mx*mx/sx2)
        sd_slope = see * s.sqrt(1./sx2)

    results=s.zeros(5)

    results[0]=slope
    results[1]=intercept
    results[2]=r
    if (len(x)>2):
        results[3]=sd_slope
        results[4]=sd_intercept

    return results







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



#y0=pylab.load('population_binned') # Now I load the initial state of the system


y0=s.zeros(len(x))
y0[0]=0.01


print 'the total population is',s.sum(y0)


N_0=s.sum(y0[:])

#I plot the initial state condition







# part dealing with the kernels

#but first the monodisperse analytical approximation

def monodisperse_approx(N_0,t,k):
	N_of_t=N_0/(1+k*N_0*t/2.)
	return N_of_t



#classical expression for the brownian kernel in the continuum regime but without the 
# slip correction and not allowing for T as a varying quantity

def Brow_ker_cont_optim(Vlist):
	kern_mat=2.*k_B*T_0/(3.*mu)*(Vlist[:,s.newaxis]**(1./df_tot)+\
	Vlist[s.newaxis,:]**(1./df_tot))**2./ \
	(Vlist[:,s.newaxis]**(1./df_tot)*Vlist[s.newaxis,:]**(1./df_tot))
	return kern_mat

def constant_kernel(Vlist):
	my_vec=Vlist/Vlist  #a vector of 1's with the right length
	kern_mat=my_vec[:,s.newaxis]*my_vec[s.newaxis,:]*8.*k_B*T_0/(3.*mu)
	return kern_mat



# def Brow_ker_cont_optim_diffusion(Vlist):
# 	#monomer volume
# 	r_mon=0.5 #monomer radius
# 	#v_mon=4./3.*s.pi*r_mon**3.
# 	#v_mono already defined as a global variable 
# 	n_mon=Vlist/v_mono #number of monomers in each aggregate

# 	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	
# 	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
# 		  (Vlist[:,s.newaxis]**(1./D_f)+Vlist[s.newaxis,:]**(1./D_f))
# 	return kern_mat


def Brow_ker_cont_optim_diffusion_adjusted(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	
	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration
	
	

	R_list[small]=a_small*n_mon[small]**(1./df_small)

	R_list[large]=a_large*n_mon[large]**(1./df_large)

	R_list[0]=0.5   #special case for the monomer radius

	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])


	return kern_mat



def Brow_ker_cont_optim_diffusion_adjusted_single_fit(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	#n_mon=Vlist/v_mono #number of monomers in each aggregate

	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration
	
# 	small=s.where(n_mon<=15.)
# 	large=s.where(n_mon>15.)

# 	a_small=0.36
# 	df_small=2.20

# 	a_large=0.24
# 	df_large=1.62

#         a_tot=0.32

# 	df_tot=1.87

	R_list=a_tot*n_mon**(1./df_tot)

#	R_list[large]=a_small*n_mon[large]**(1./df_large)

	R_list[0]=0.5   #special case for the monomer radius

	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])


	return kern_mat





def Brow_ker_cont_optim_diffusion_adjusted_and_monomer(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	#n_mon=Vlist/v_mono #number of monomers in each aggregate

	print "n_mon is, ", n_mon
	
	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration

	#threshold=10.
	
# 	small=s.where(n_mon<=threshold)
# 	large=s.where(n_mon>threshold)

# 	a_small=0.385
# 	df_small=2.417

# 	a_large=0.245
# 	df_large=1.63


	R_list[small]=a_small*n_mon[small]**(1./df_small)

	R_list[large]=a_large*n_mon[large]**(1./df_large)

#	R_list[0]=0.5   #special case for the monomer radius

	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])


	return kern_mat




def Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	#n_mon=Vlist/v_mono #number of monomers in each aggregate

	#print "n_mon is, ", n_mon
	
	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration

	#threshold=15.
	
# 	small=s.where(n_mon<=threshold)
# 	large=s.where(n_mon>threshold)

# 	a_small=0.36
# 	df_small=2.19

# 	a_large=0.241
# 	df_large=1.62


	R_list[small]=a_small*n_mon[small]**(1./df_small)

	R_list[large]=a_large*n_mon[large]**(1./df_large)

#	R_list[0]=0.5   #special case for the monomer radius

#	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

        #m=Vlist #since rho = 1.

	c=(8.*k_B*T_0/(s.pi*m))**0.5
	#print 'c is', c
	l=8.*Diff/(s.pi*c)
	#print 'l is', l
	diam_seq=2.*R_list
	
	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
	
	beta_fuchs=(\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
	)**(-1.)
	
	
	## now I have all the bits for the kernel matrix
#	kern_mat=Brow_ker_cont_optim(Vlist)*beta


	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])*beta_fuchs


	return kern_mat



def Brow_ker_cont_optim_diffusion_adjusted_and_monomer_single_fit_beta_fuchs(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	#n_mon=Vlist/v_mono #number of monomers in each aggregate

	#print "n_mon is, ", n_mon
	
	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration

# 	threshold=15.
	
# 	small=s.where(n_mon<=threshold)
# 	large=s.where(n_mon>threshold)

# 	a_small=0.36
# 	df_small=2.19

# 	a_large=0.241
# 	df_large=1.62


# 	R_list[small]=a_small*n_mon[small]**(1./df_small)

# 	R_list[large]=a_small*n_mon[large]**(1./df_large)



	R_list=a_tot*n_mon**(1./df_tot)



#	R_list[0]=0.5   #special case for the monomer radius

#	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

        #m=Vlist #since rho = 1.

	c=(8.*k_B*T_0/(s.pi*m))**0.5
	#print 'c is', c
	l=8.*Diff/(s.pi*c)
	#print 'l is', l
	diam_seq=2.*R_list
	
	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
	
	beta_fuchs=(\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
	)**(-1.)
	
	
	## now I have all the bits for the kernel matrix
#	kern_mat=Brow_ker_cont_optim(Vlist)*beta


	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])*beta_fuchs


	return kern_mat





def Brow_ker_cont_optim_standard_diffusion_adjusted_and_monomer_beta_fuchs(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	#n_mon=Vlist/v_mono #number of monomers in each aggregate

	#print "n_mon is, ", n_mon
	

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration



	#threshold=15.
	
# 	small=s.where(n_mon<=threshold)
# 	large=s.where(n_mon>threshold)

# 	a_small=0.36
# 	df_small=2.19

# 	a_large=0.241
# 	df_large=1.62


	R_list[small]=a_small*n_mon[small]**(1./df_small)

	R_list[large]=a_large*n_mon[large]**(1./df_large)

	Diff=k_B*T_0/(6.*s.pi*mu*R_list) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.


#	R_list[0]=0.5   #special case for the monomer radius

#	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

        #m=Vlist #since rho = 1.

	c=(8.*k_B*T_0/(s.pi*m))**0.5
	#print 'c is', c
	l=8.*Diff/(s.pi*c)
	#print 'l is', l
	diam_seq=2.*R_list
	
	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
	
	beta_fuchs=(\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
	)**(-1.)

	print "beta_fuchs is, ", beta_fuchs
	
	## now I have all the bits for the kernel matrix
#	kern_mat=Brow_ker_cont_optim(Vlist)*beta


	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])*beta_fuchs


	return kern_mat



def Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_single_value(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	#n_mon=Vlist/v_mono #number of monomers in each aggregate

	#print "n_mon is, ", n_mon
	
	

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration



# 	threshold=15.
	
# 	small=s.where(n_mon<=threshold)
# 	large=s.where(n_mon>threshold)

# 	a_small=0.36
# 	df_small=2.19

# 	a_large=0.241
# 	df_large=1.62


# 	R_list[small]=a_small*n_mon[small]**(1./df_small)

# 	R_list[large]=a_small*n_mon[large]**(1./df_large)



#	R_list=a_tot*n_mon**(1./df_tot)

	R_list[small]=a_small*n_mon[small]**(1./df_small)

	R_list[large]=a_large*n_mon[large]**(1./df_large)



	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.

        print "the diffusion coefficient is, ", Diff


#	R_list[0]=0.5   #special case for the monomer radius

#	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

#        m=Vlist #since rho = 1.

# 	c=(8.*k_B*T_0/(s.pi*m))**0.5
# 	#print 'c is', c
# 	l=8.*Diff/(s.pi*c)
# 	#print 'l is', l
# 	diam_seq=2.*R_list
	
# 	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
	
# 	beta_fuchs=(\
# 	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
# 	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
# 	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
# 	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
# 	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
# 	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
# 	)**(-1.)
	
	
	## now I have all the bits for the kernel matrix
#	kern_mat=Brow_ker_cont_optim(Vlist)*beta


        #beta_fuchs=beta_fuchs_standalone(m_mon)

        #beta_fuchs=beta_fuchs[0,0]  #I am correcting for monomers only

	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])*beta_00


	return kern_mat



def Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_special_rules(Vlist):
	#monomer volume
	#r_mon=0.5 #monomer radius
	#v_mon=4./3.*s.pi*r_mon**3.
	#v_mono already defined as a global variable 
	#n_mon=Vlist/v_mono #number of monomers in each aggregate

	#print "n_mon is, ", n_mon
	
	Diff=k_B*T_0/(n_mon*m_mon*beta) #vector with the cluster diffusion coefficients
	#which just depends on the cluster size.

	R_list=s.zeros(len(Vlist)) #I initialize the array which will contain the
	#radia of gyration

	#threshold=15.
	
# 	small=s.where(n_mon<=threshold)
# 	large=s.where(n_mon>threshold)

# 	a_small=0.36
# 	df_small=2.19

# 	a_large=0.241
# 	df_large=1.62


	R_list[small]=a_small*n_mon[small]**(1./df_small)

	R_list[large]=a_large*n_mon[large]**(1./df_large)

        #Now I actually introduce some special rules for the radius of gyration of the smallest clusters

        sel=s.where(n_mon == 1.)

        R_list[sel]=3.87e-1


        sel=s.where(n_mon == 2.)

        R_list[sel]=6.39e-1

        sel=s.where(n_mon == 3.)

        R_list[sel]=7.18e-1

        sel=s.where(n_mon == 4.)

        R_list[sel]=7.48e-1


        print "len(R_list is, )", len(R_list)
        print "R_list is, ", R_list


        

#	R_list[0]=0.5   #special case for the monomer radius

#	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

        #m=Vlist #since rho = 1.

	c=(8.*k_B*T_0/(s.pi*m))**0.5
	#print 'c is', c
	l=8.*Diff/(s.pi*c)
	#print 'l is', l
	diam_seq=2.*R_list
	
	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
	
	beta_fuchs=(\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
	)**(-1.)
	
	
	## now I have all the bits for the kernel matrix
#	kern_mat=Brow_ker_cont_optim(Vlist)*beta


	kern_mat=4.*s.pi*(Diff[:,s.newaxis]+Diff[s.newaxis,:])* \
		  (R_list[:,s.newaxis]+R_list[s.newaxis,:])*beta_fuchs


	return kern_mat



    
def beta_fuchs_standalone(Vlist):

        #m=Vlist #since rho = 1.


        n_mon=1  #I am calculating the correct diffusion coefficient for a single monomer

        print "n_mon is,", n_mon

	Diff=k_B*T_0/(n_mon*m_mon*beta)

	c=(8.*k_B*T_0/(s.pi*m))**0.5
	#print 'c is', c
	l=8.*Diff/(s.pi*c)
	#print 'l is', l
        R_list=0.5 #radius of a monomer
	diam_seq=2.*R_list
	
	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq

	beta_fuchs=(\
	(diam_seq+diam_seq)   \
	/(diam_seq+diam_seq\
	+2.*(g**2.+g**2.)**0.5)\
	+ 8.*(Diff+Diff)/\
	((c**2.+c**2.)**0.5*\
	(diam_seq+diam_seq))\
	)**(-1.)

        return beta_fuchs

        
	




def coupling_optim_garrick(y,t):
	#now the kernel is a global variable like lensum
	creation=s.zeros(lensum)
	destruction=s.zeros(lensum)
	#now I try to rewrite this in a more optimized way
	destruction = -s.dot(s.transpose(kernel),y)*y #much more concise way to express\
	#the destruction of k-mers 
	
	for k in xrange(lensum):
		kyn = (kernel*f_garrick[:,:,k])*y[:,s.newaxis]*y[s.newaxis,:]
		creation[k] = s.sum(kyn)
	creation=0.5*creation
	out=creation+destruction
	return out



#Now I work with the function for espressing smoluchowski equation when a uniform grid is used

def coupling_optim(y,t):
	creation=s.zeros(lensum)
	destruction=s.zeros(lensum)
	#now I try to rewrite this in a more optimized way
	destruction = -s.dot(s.transpose(kernel),y)*y #much more concise way to express\
	#the destruction of k-mers 
	kyn = kernel*y[:,s.newaxis]*y[s.newaxis,:]
	for k in xrange(lensum):
		creation[k] = s.sum(kyn[s.arange(k),k-s.arange(k)-1])
	creation=0.5*creation
	out=creation+destruction
	return out






#Now I go for the optimal optimization of the chi_{i,j,k} coefficients used by Garrick for
# dealing with a non-uniform grid. 

def mycount_garrick(V):
	f=s.zeros((lensum, lensum, lensum))
	#f_sat=zeros(lensum)
	Vsum=V[:,s.newaxis]+V[s.newaxis,:] # matrix with the sum of the volumes in the bins
	for k in xrange(1,(lensum-1)):
		f[:,:,k]=s.where((Vsum<=V[k+1]) & (Vsum>=V[k]), (V[k+1]-Vsum)/(V[k+1]-V[k]),\
		f[:,:,k] )
		f[:,:,k]=s.where((Vsum<=V[k]) & (Vsum>=V[k-1]),(Vsum-V[k-1])/(V[k]-V[k-1]),\
		f[:,:,k])
	
	return f	

#####################################################################################
######################################################################################
######################################################################################




beta_00 = beta_fuchs_standalone(m_mon)

print "beta_00 is, ", beta_00



if (choice_kernel == 0):
    kernel=Brow_ker_cont_optim(x)
elif (choice_kernel == 1):
    kernel=constant_kernel(x)
elif (choice_kernel == 2):
    Brow_ker_cont_optim_diffusion_adjusted(x)
elif (choice_kernel == 3):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_single_fit(x)
elif (choice_kernel == 4):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer(x)
elif (choice_kernel == 5):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs(x)
elif (choice_kernel == 6):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_single_fit_beta_fuchs(x)
elif (choice_kernel == 7):
    kernel=Brow_ker_cont_optim_standard_diffusion_adjusted_and_monomer_beta_fuchs(x)
    
elif (choice_kernel == 8):
    kernel=Brow_ker_cont_optim_standard_diffusion_adjusted_and_monomer_single_fit_beta_fuchs(x)
    

elif (choice_kernel == 9):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_single_value(x)

elif (choice_kernel ==10):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_special_rules(x)


f_garrick=mycount_garrick(x)

print "ok down to here"





y = si.odeint(coupling_optim_garrick, y0, \
		   t,printmessg=1,rtol=1e-10,atol=1e-10)
	


p.save("evolution_of_the_population.dat",y)

concentration=s.sum(y,axis=1)


count_clus=n_part/rho

concentration=concentration*count_clus

print "the number of monomer is, ", concentration

p.save("number_clusters_smoluchowski.dat", concentration)

p.save("time_smoluchowski.dat",t)


mass_in_time=(y*x).sum(axis=1)

p.save("mass_conservation.dat", mass_in_time)

#Now I discuss the mear R and the mean number of monomers in a cluster

k_aver=n_part/concentration

p.save("average_number_monomers_in_cluster.dat", k_aver)




fig = p.figure()
axes = fig.gca()

#axes.plot(s.log10(n_comp_large),s.log10(r_gyr_ari_large), "bo")
axes.plot(t,k_aver, "ro",label="mean number of monomers in a cluster")
#axes.plot(t,number_mono,label="monodisperse approx",linewidth=2.)
# axes.loglog(time_clu,n_clus_large, "rx",label="$N_>$")
# axes.loglog(time_clu[choice_time_clu],(10.**my_conc),"k",\
#             label=r"Fit of $N_\infty $",linewidth=2.)
p.xlabel('Time')
p.ylabel('Number of monomers in cluster')
p.title("Evolution Number of monomers in a  clusters")
p.grid(True)
cluster_name="number_of_monomers_in_a_cluster.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()




k_select_small=s.where(k_aver<=threshold)
k_select_large=s.where(k_aver>threshold)




R_aver=s.zeros(len(k_aver))

R_aver[k_select_small]=a_small*k_aver[k_select_small]**(1./df_small)

R_aver[k_select_large]=a_small*k_aver[k_select_large]**(1./df_large)


p.save("average_radius_of_gyration_smoluchowski_wrong.dat", R_aver)



#R_aver2=s.zeros(len(k_aver))

# k_vec=x/v_mono


# print "k_vec is, ", k_vec

#just a test; I am recalculating the same object twice.

k_aver2=(y*n_mon).sum(axis=1)/y.sum(axis=1)


p.save("average_number_monomers_in_cluster.dat", k_aver2)







R_vec=s.zeros(len(n_mon))

if (choice_slope == 2):

    R_vec[small]=a_small*n_mon[small]**(1./df_small)

    R_vec[large]=a_large*n_mon[large]**(1./df_large)

elif (choice_slope == 1):
    R_vec=a_tot*n_mon**(1./df_tot)


R_aver2=(y*R_vec).sum(axis=1)/y.sum(axis=1)


p.save("average_radius_of_gyration_smoluchowski2.dat", R_aver2)


fig = p.figure()
axes = fig.gca()

#axes.plot(s.log10(n_comp_large),s.log10(r_gyr_ari_large), "bo")
axes.plot(t,R_aver2, "ro",label="mean number cluster R_g")
#axes.plot(t,number_mono,label="monodisperse approx",linewidth=2.)
# axes.loglog(time_clu,n_clus_large, "rx",label="$N_>$")
# axes.loglog(time_clu[choice_time_clu],(10.**my_conc),"k",\
#             label=r"Fit of $N_\infty $",linewidth=2.)
p.xlabel('Time')
p.ylabel('Mean R_g in cluster 2')
p.title("Evolution mean radius of gyration 2")
p.grid(True)
cluster_name="mean_radius_of_gyration_in_a_cluster2.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()









fig = p.figure()
axes = fig.gca()

#axes.plot(s.log10(n_comp_large),s.log10(r_gyr_ari_large), "bo")
axes.plot(t,k_aver2, "ro",label="mean number cluster R_g")
#axes.plot(t,number_mono,label="monodisperse approx",linewidth=2.)
# axes.loglog(time_clu,n_clus_large, "rx",label="$N_>$")
# axes.loglog(time_clu[choice_time_clu],(10.**my_conc),"k",\
#             label=r"Fit of $N_\infty $",linewidth=2.)
p.xlabel('Time')
p.ylabel('Mean k in cluster')
p.title("Evolution mean k in a cluster")
p.grid(True)
cluster_name="number_of_monomers_in_a_cluster2.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()








fig = p.figure()
axes = fig.gca()

#axes.plot(s.log10(n_comp_large),s.log10(r_gyr_ari_large), "bo")
axes.plot(t,R_aver, "ro",label="mean number cluster R_g")
#axes.plot(t,number_mono,label="monodisperse approx",linewidth=2.)
# axes.loglog(time_clu,n_clus_large, "rx",label="$N_>$")
# axes.loglog(time_clu[choice_time_clu],(10.**my_conc),"k",\
#             label=r"Fit of $N_\infty $",linewidth=2.)
p.xlabel('Time')
p.ylabel('Mean R_g in cluster')
p.title("Evolution mean radius of gyration")
p.grid(True)
cluster_name="mean_radius_of_gyration_in_a_cluster_wrong.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()





#now the monodisperse approx

kernel_const=8.*k_B*T_0/(3.*mu)

number_mono=monodisperse_approx(0.01,t,kernel_const)

number_mono=number_mono*count_clus

p.save("number_mono_analytic.dat", number_mono)



fig = p.figure()
axes = fig.gca()

#axes.plot(s.log10(n_comp_large),s.log10(r_gyr_ari_large), "bo")
axes.plot(t,concentration, "ro",label="concentration numerics")
axes.plot(t,number_mono,label="monodisperse approx",linewidth=2.)
# axes.loglog(time_clu,n_clus_large, "rx",label="$N_>$")
# axes.loglog(time_clu[choice_time_clu],(10.**my_conc),"k",\
#             label=r"Fit of $N_\infty $",linewidth=2.)
p.xlabel('Time')
p.ylabel('Number of clusters')
p.title("Evolution Number of clusters")
p.grid(True)
cluster_name="test_constant_kernel_and_monodisperse.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()



time_sel=s.where(t>=7000.)


lin_fit=fitting_stat(s.log10(t[time_sel]),s.log10(concentration[time_sel]))

print "the slope is, ", lin_fit[0]




print "Now ahead with another integration, this time on a linear bin structure"

#I have to re-define a number of important quantities


n_mon=s.linspace(1., 2500., 2500)  #new volume grid


x=n_mon*v_mono  #new volume grid



print "n_mon is, ", n_mon

m=n_mon #since each monomer has mass 1 in my units

small=s.where(n_mon<=threshold)
large=s.where(n_mon>threshold)


lensum=len(x) # number of particle bins

#and also the initial state is different since it is defined on a new grid

y0=s.zeros(len(x))
y0[0]=0.01


#ricalculate kernel matrix on new volume grid



if (choice_kernel == 0):
    kernel=Brow_ker_cont_optim(x)
elif (choice_kernel == 1):
    kernel=constant_kernel(x)
elif (choice_kernel == 2):
    Brow_ker_cont_optim_diffusion_adjusted(x)
elif (choice_kernel == 3):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_single_fit(x)
elif (choice_kernel == 4):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer(x)
elif (choice_kernel == 5):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs(x)
elif (choice_kernel == 6):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_single_fit_beta_fuchs(x)
elif (choice_kernel == 7):
    kernel=Brow_ker_cont_optim_standard_diffusion_adjusted_and_monomer_beta_fuchs(x)
    
elif (choice_kernel == 8):
    kernel=Brow_ker_cont_optim_standard_diffusion_adjusted_and_monomer_single_fit_beta_fuchs(x)
    

elif (choice_kernel == 9):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_single_value(x)

elif (choice_kernel ==10):
    kernel=Brow_ker_cont_optim_diffusion_adjusted_and_monomer_beta_fuchs_special_rules(x)


y_new = si.odeint(coupling_optim, y0, \
		   t,printmessg=1,rtol=1e-10,atol=1e-10)



p.save("evolution_of_the_population_linear_bins.dat",y_new)

concentration=s.sum(y_new,axis=1)

count_clus=5000./0.01

concentration=concentration*count_clus

p.save("number_clusters_smoluchowski_linear_bins.dat", concentration)





# now a test on mass conservation

mass_new_in_time=(y_new*x).sum(axis=1)

p.save("mass_new_conservation_linear_bins.dat", mass_new_in_time)


# delta_vol=x[30]-x[29]




# print "the mass evolution in time is, ", mass_new_in_time






print 'So far so good'