myprojects-hg-reborn

Fork

#! /usr/bin/env python

import scipy as s
import numpy as n
import pylab as p
#from rpy import r
import distance_calc as d_calc
import igraph as ig
import hydrodynamic_radius_calc as h
import calc_rad as c_rad

from scipy import stats #I need this module for the linear fit


box_size=10000.  #here the box size does not depend on the total particle number
#but on the total number of particles in a box
print "the box size is, ", box_size

#NB: now the box size is fixed a priori and is no longer a derivate object

n_part=5000
#density=0.03

size_cluster=50 #number of monomers in each cluster.

types=n_part/size_cluster #different kinds of particles i.e. total number of clusters when coagulation
# is done

print "types is, ", types


min_clus_stat=250  # this defines the minimum number of clusters which need to have formed before I start taking some statistics


# x_cm=s.zeros(types)   
# y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
# z_cm=s.zeros(types)

# v_x_cm=s.zeros(types)
# v_y_cm=s.zeros(types)  #velocities of the centres of mass of each cluster
# v_z_cm=s.zeros(types)


#Again, I prefer to give the number of particles and the density as input parameters
# and work out the box size accordingly


ini_config=0
fin_config=1500 #for large data post-processing I need to declare an initial and final

               #configuration I want to read and post-process


by=10 #this tells how many configurations there are in the file I am reading

figure=0 #whether I sould print many figures or not


n_config=(fin_config-ini_config)/by # total number of configurations, which are numbered from 0 to n_config-1
my_selection=s.arange(ini_config,fin_config,by)


threshold=1.06 #threshold to consider to particles as directly connected




time=p.load("time.dat")

time=time[my_selection] #I adapt the time to the initial and final configuration

p.save("time_red.dat",time)

# x_mean=s.zeros(len(time))     # mean position (average on all the clusters) of the cluster CM
# y_mean=s.zeros(len(time))
# z_mean=s.zeros(len(time))

# #y_sq_mean=s.zeros(len(time))
# #z_sq_mean=s.zeros(len(time))
# #z_sq_mean=s.zeros(len(time))


# delta_x_sq=s.zeros(len(time))
# delta_y_sq=s.zeros(len(time))
# delta_z_sq=s.zeros(len(time))



# # v_x_sq_mean=s.zeros(len(time))
# # v_y_sq_mean=s.zeros(len(time))
# # v_z_sq_mean=s.zeros(len(time))

# delta_v_x_sq=s.zeros(len(time))
# delta_v_y_sq=s.zeros(len(time))
# delta_v_z_sq=s.zeros(len(time))



#dist_mat=s.zeros((size_cluster,size_cluster)) # I need to find out when the particles of each kind
#coagulate into a single cluster...I am not interested in looking at at all the particles in the system now



def R_cm_distr(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
    #and velocity

    count_left=s.arange(0,n_part,size_cluster)
    count_right=count_left+size_cluster         


    
#     x_cm=s.zeros(types)   
#     y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
#     z_cm=s.zeros(types)

    
#     #print "OK here"
#     x_pos_temp=test_arr[:,0]
#     y_pos_temp=test_arr[:,1]
#     z_pos_temp=test_arr[:,2]

    R_cm_arr=s.zeros((types,3)) #this array will contain the distribution of hydrodynamic radia
    #calculated for a certain snapshot of the system


    for i in xrange(types):
        r_pos=test_arr[count_left[i]:count_right[i],:]
#         x_pos=x_pos_temp[count_left[i]:count_right[i]]
#         y_pos=y_pos_temp[count_left[i]:count_right[i]]
#         z_pos=z_pos_temp[count_left[i]:count_right[i]]

     

#        print "s.shape(r_pos) is, ", s.shape(r_pos)
        
        R_cm_arr[i,:]=r_pos.mean(axis=0)


    return R_cm_arr






def R_h_distr(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
    #and velocity

    count_left=s.arange(0,n_part,size_cluster)
    count_right=count_left+size_cluster         


    
    x_cm=s.zeros(types)   
    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
    z_cm=s.zeros(types)

    
    #print "OK here"
    x_pos_temp=test_arr[:,0]
    y_pos_temp=test_arr[:,1]
    z_pos_temp=test_arr[:,2]

    R_hydro_arr=s.zeros(types) #this array will contain the distribution of hydrodynamic radia
    #calculated for a certain snapshot of the system


    for i in xrange(types):
        x_pos=x_pos_temp[count_left[i]:count_right[i]]
        y_pos=y_pos_temp[count_left[i]:count_right[i]]
        z_pos=z_pos_temp[count_left[i]:count_right[i]]

        #print "x_pos is, ", x_pos
        
        R_hydro_arr[i]=h.hydro_radius_calc(x_pos,y_pos,z_pos,box_size)
   #     print "R_hydro_arr[i] is, ", R_hydro_arr[i]
        
    R_h=R_hydro_arr.mean()


    return R_h





def R_gyr_distr(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
    #and velocity

    count_left=s.arange(0,n_part,size_cluster)
    count_right=count_left+size_cluster         


    
    x_cm=s.zeros(types)   
    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
    z_cm=s.zeros(types)

    
    #print "OK here"
    x_pos_temp=test_arr[:,0]
    y_pos_temp=test_arr[:,1]
    z_pos_temp=test_arr[:,2]

    R_gyr_arr=s.zeros(types) #this array will contain the distribution of hydrodynamic radia
    #calculated for a certain snapshot of the system


    for i in xrange(types):
        x_pos=x_pos_temp[count_left[i]:count_right[i]]
        y_pos=y_pos_temp[count_left[i]:count_right[i]]
        z_pos=z_pos_temp[count_left[i]:count_right[i]]

        #print "x_pos is, ", x_pos
        
        R_gyr_arr[i]=c_rad.rad_gyr(x_pos,y_pos,z_pos,box_size)
   #     print "R_hydro_arr[i] is, ", R_hydro_arr[i]
        
    R_gyr=R_gyr_arr.mean()


    return R_gyr




def calc_time_stat(test_arr,types,n_config,n_part,size_cluster,box_size):  #I can use this function both for the statistics on position
    #and velocity

    count_left=s.arange(0,n_part,size_cluster)
    count_right=count_left+size_cluster         


    
    x_cm=s.zeros(types)   
    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
    z_cm=s.zeros(types)


    #fold everything inside the box

    test_arr=s.remainder(test_arr,box_size)
    
    #print "OK here"
    x_pos=test_arr[:,0]
    y_pos=test_arr[:,1]
    z_pos=test_arr[:,2]


    results=s.zeros(7)


    for j in xrange(types):
        x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
        y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
        z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])





    x_mean=x_cm.mean()
    y_mean=y_cm.mean()
    z_mean=z_cm.mean()

    delta_x_sq=x_cm.var()
    delta_y_sq=y_cm.var()
    delta_z_sq=z_cm.var()

    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq


    results[0]=x_mean
    results[1]=y_mean
    results[2]=z_mean

    results[3]=delta_x_sq
    results[4]=delta_y_sq
    results[5]=delta_z_sq
    results[6]=delta_r_sq


    return results

def calc_time_stat_vel(test_arr,types,n_config,n_part,size_cluster,box_size):  #I can use this function both for the statistics on position
    #and velocity

    count_left=s.arange(0,n_part,size_cluster)
    count_right=count_left+size_cluster         


    
    x_cm=s.zeros(types)   
    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
    z_cm=s.zeros(types)

#here I have the only difference wrt the position fluctuations: I do not have to fold the velocity!!!
 #    #fold everything inside the box

#     test_arr=s.remainder(test_arr,box_size)
    
    #print "OK here"
    x_pos=test_arr[:,0]
    y_pos=test_arr[:,1]
    z_pos=test_arr[:,2]


    results=s.zeros(7)


    for j in xrange(types):
        x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
        y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
        z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])





    x_mean=x_cm.mean()
    y_mean=y_cm.mean()
    z_mean=z_cm.mean()

    delta_x_sq=x_cm.var()
    delta_y_sq=y_cm.var()
    delta_z_sq=z_cm.var()

    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq


    results[0]=x_mean
    results[1]=y_mean
    results[2]=z_mean

    results[3]=delta_x_sq
    results[4]=delta_y_sq
    results[5]=delta_z_sq
    results[6]=delta_r_sq


    return results








def calc_time_stat_clusters_together(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
    #and velocity

    count_left=s.arange(0,n_part,size_cluster)
    count_right=count_left+size_cluster         


    
    x_cm=s.zeros(types)   
    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
    z_cm=s.zeros(types)

    
    #print "OK here"
    x_pos=test_arr[:,0]
    y_pos=test_arr[:,1]
    z_pos=test_arr[:,2]


    results=s.zeros(7)







    x_mean=x_pos.mean()
    y_mean=y_pos.mean()
    z_mean=z_pos.mean()

    delta_x_sq=x_pos.var()
    delta_y_sq=y_pos.var()
    delta_z_sq=z_pos.var()

    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq


    results[0]=x_mean
    results[1]=y_mean
    results[2]=z_mean

    results[3]=delta_x_sq
    results[4]=delta_y_sq
    results[5]=delta_z_sq
    results[6]=delta_r_sq


    return results












def calc_time_stat_selection(test_arr,types,n_config,n_part,size_cluster,clust_choice):
#I can use this function both for the statistics on position
    #and velocity


    
    x_cm=s.zeros(len(s.ravel(clust_choice)))   
    y_cm=s.zeros(len(s.ravel(clust_choice)))    #positions of the centres of mass of each cluster
    z_cm=s.zeros(len(s.ravel(clust_choice)))


   # print "len(clust_choice) is, ", len(clust_choice)
    count_left=s.arange(0,(len(s.ravel(clust_choice))*size_cluster),size_cluster)
    count_right=count_left+size_cluster         



    
    #print "OK here"
    x_pos_temp=test_arr[:,0]
    y_pos_temp=test_arr[:,1]
    z_pos_temp=test_arr[:,2]

    #print "len(x_pos_temp) is, ", len(x_pos_temp)

    x_pos=s.zeros(len(s.ravel(clust_choice))*size_cluster)
    y_pos=s.zeros(len(s.ravel(clust_choice))*size_cluster)
    z_pos=s.zeros(len(s.ravel(clust_choice))*size_cluster)

    #print "len(x_pos) is,", len(x_pos)

    for l in xrange(len(s.ravel(clust_choice))):
     #   print "count_left and count_right are, ", count_left[l],count_right[l]

        #print "clust_choice[l] is, ", clust_choice[l]
        x_pos[count_left[l]:count_right[l]]=\
        x_pos_temp[(s.ravel(clust_choice)[l]*size_cluster):(s.ravel(clust_choice)[l]*size_cluster+size_cluster)]
        y_pos[count_left[l]:count_right[l]]=\
                 y_pos_temp[(s.ravel(clust_choice)[l]*size_cluster):(s.ravel(clust_choice)[l]*size_cluster+size_cluster)]

        z_pos[count_left[l]:count_right[l]]=\
                 z_pos_temp[(s.ravel(clust_choice)[l]*size_cluster):(s.ravel(clust_choice)[l]*size_cluster+size_cluster)]





    
    #Then all the rest can stay the same

    results=s.zeros(7)



    for j in xrange(len(s.ravel(clust_choice))): #I iterate on the existing fully-coagulated clusters
        #print "j is, ", j
        x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
        y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
        z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])





    x_mean=x_cm.mean()
    y_mean=y_cm.mean()
    z_mean=z_cm.mean()

    delta_x_sq=x_cm.var()
    delta_y_sq=y_cm.var()
    delta_z_sq=z_cm.var()

    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq


    results[0]=x_mean
    results[1]=y_mean
    results[2]=z_mean

    results[3]=delta_x_sq
    results[4]=delta_y_sq
    results[5]=delta_z_sq
    results[6]=delta_r_sq


    return results







def coagulation_clusters(test_arr,types,n_config,box_size,size_cluster,threshold):

    count_left=s.arange(0,n_part,size_cluster)
    count_right=count_left+size_cluster         
    
    coag_arr=s.zeros(types)

    
#     x_cm=s.zeros(types)   
#     y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
#     z_cm=s.zeros(types)

    
    #print "OK here"
    x_pos=test_arr[:,0]
    y_pos=test_arr[:,1]
    z_pos=test_arr[:,2]


    results=s.zeros(7)


    for j in xrange(types):
        x_pos_clu=x_pos[count_left[j]:count_right[j]]
        y_pos_clu=y_pos[count_left[j]:count_right[j]]
        z_pos_clu=z_pos[count_left[j]:count_right[j]]

        dist_mat=d_calc.distance_calc(x_pos_clu,y_pos_clu,z_pos_clu, box_size)
        cluster_obj=ig.Graph.Adjacency((dist_mat <= threshold).tolist(),\
					       ig.ADJ_UNDIRECTED)

        cluster_obj.simplify()
        clustering=cluster_obj.clusters()
        
        #n_clusters[i]=p.load("nc_temp.dat")
#        print "clustering is, ", clustering

        coag_arr[j]=len(clustering)

   # print "coag_arr is, ", coag_arr
    return coag_arr


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







# #print "Ok here"

# for i in xrange(0,n_config):

#     read_pos="read_pos_%1d"%my_selection[i]
#     read_vel="read_vel_%1d"%my_selection[i]

#     #print "OK here 2"
#     #cluster_name="cluster_dist%05d"%my_selection[i]
#     test_arr=p.load(read_pos)
    
#     #print "OK here"
#     x_pos=test_arr[:,0]
#     y_pos=test_arr[:,1]
#     z_pos=test_arr[:,2]

#     for j in xrange(types):
#         x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
#         y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
#         z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])

#         v_x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
#         v_y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
#         v_z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])



#     x_mean[i]=x_cm.mean()
#     y_mean[i]=y_cm.mean()
#     z_mean[i]=z_cm.mean()

#     delta_x_sq[i]=x_cm.var()
#     delta_y_sq[i]=y_cm.var()
#     delta_z_sq[i]=z_cm.var()



# delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq


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



# axes.plot(time,x_mean,"ro",\
#           label="<x>",linewidth=2.)
# axes.plot(time,y_mean,"k+",\
#           label="<y>")
# axes.plot(time,z_mean,"mx",\
#           label="<z> ",linewidth=2.)
# #p.ylim(1.4,2.6)
# #axes.legend()
# p.xlabel('Time')
# p.ylabel('Mean position')
# p.savefig("mean_x_y_z.pdf")
 
# p.clf()    

        
   

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



# axes.plot(time,delta_r_sq,"ro",\
#           label="<x>",linewidth=2.)
# # axes.plot(time,y_mean,"k+",\
# #           label="<y>")
# # axes.plot(time,z_mean,"mx",\
# #           label="<z> ",linewidth=2.)
# #p.ylim(1.4,2.6)
# #axes.legend()
# p.xlabel('Time')
# p.ylabel('delta position squared')
# p.savefig("delta_r_squared.pdf")
 
# p.clf()    


    
#     #print "the length of x_pos is, ", len(x_pos)

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

#initialization of some arrays I will have to take in & out the loop


#pos_results=s.zeros(4) #just to initialize it
#pos_results_selection=s.zeros(4) #just to initialize it
#pos_statistics_selection=s.zeros(2)  # I just need to initialize it!
time_selection=s.zeros(2)




flag_ini=0 # I will need this flag to check is some other array has been initialized

R_h_arr=s.zeros(n_config)

R_gyr_arr=s.zeros(n_config)


for i in xrange(0,n_config):

    read_pos="read_pos_%1d"%my_selection[i]
    read_vel="read_vel_%1d"%my_selection[i]

    #print "OK here 2"
    #cluster_name="cluster_dist%05d"%my_selection[i]
    pos_arr=p.load(read_pos)
    vel_arr=p.load(read_vel)
        
    pos_statistics=calc_time_stat(pos_arr,types,n_config,n_part,size_cluster, box_size)        

#    pos_statistics_clusters_together=calc_time_stat_clusters_together(pos_arr,types,n_config,n_part,size_cluster)        


    vel_statistics=calc_time_stat_vel(vel_arr,types,n_config,n_part,size_cluster,box_size)        
    clust_stat=coagulation_clusters(pos_arr,types,n_config,box_size,size_cluster,threshold)
    #print "clust_stat is, ", clust_stat
    clust_choice=s.where(clust_stat[:]==1.) #i.e. I select the clusters which have already finished
    # coagulating

    R_h_arr[i]=R_h_distr(pos_arr,types,n_config,n_part,size_cluster)
    R_gyr_arr[i]=R_gyr_distr(pos_arr,types,n_config,n_part,size_cluster)
    R_cm=R_cm_distr(pos_arr,types,n_config,n_part,size_cluster)

    if (i==0):
        R_cm_arr=R_cm
    elif (i>0):
        R_cm_arr=s.hstack((R_cm_arr,R_cm))

    print "R_gyr_arr[i] is, ", R_gyr_arr[i]
    print "R_h_arr[i] is, ", R_h_arr[i]
    
    #print "len(clust_choice) is, ", len(s.ravel(clust_choice))
    #print "clust_choice is, ", clust_choice

    if (len(s.ravel(clust_choice))>=min_clus_stat):
        pos_statistics_selection=\
             calc_time_stat_selection(pos_arr,types,n_config,n_part,size_cluster,clust_choice)
     #   print "pos_statistics_selection is, ", pos_statistics_selection

    



    if (i==0):
        
        pos_results=pos_statistics
        #R_h=R_h_temp

#        pos_results_clusters_together=pos_statistics_clusters_together


#         print "pos_results is, ", pos_results
#         print "the len of pos_results is, ", len(pos_results)
        vel_results=vel_statistics
        clust_results=clust_stat
        #The following if condition is academic since I'll never have the right number of
        #clusters at the beginning
        if (len(s.ravel(clust_choice))>=min_clus_stat):
            pos_results_selection=pos_statistics_selection
            time_selection=[time[i]]
            
        
            
    elif(i!=0):
#         print "now pos_results is, ", pos_results
#         print "now the shape of pos_results is, ", s.shape(pos_results)
#         print "now pos_statistics is, ", pos_statistics

        pos_results=s.concatenate((pos_results,pos_statistics))
        #R_h=s.concatenate((R_h,R_h_temp))
#        pos_results_clusters_together=s.concatenate((pos_results_clusters_together,pos_statistics_clusters_together))



#         print "OK here"
        vel_results=s.concatenate((vel_results,vel_statistics))
        clust_results=s.concatenate((clust_results,clust_stat))


        if ((len(s.ravel(clust_choice))>=min_clus_stat) and (flag_ini==0)):
            pos_results_selection=pos_statistics_selection
            time_selection=[time[i]]
            flag_ini=1 #to be sure I am doing this only the fist time, then I can concatenate


        elif ((len(s.ravel(clust_choice))>=min_clus_stat) and (flag_ini!=0)):
            pos_results_selection=s.concatenate((pos_results_selection,pos_statistics_selection))
            time_selection=s.concatenate((time_selection,[time[i]] ))



pos_results.shape=(-1,7)
vel_results.shape=(-1,7)  #where -1 means "do whatever is needed to reshape it right"



#pos_results_clusters_together.shape=(-1,7)





clust_results.shape=(-1,types)


#pos_results_selection.shape=(-1,7)


p.save("position_cluster_centres_of_mass.dat", pos_results)
p.save("velocity_cluster_centres_of_mass.dat", vel_results)
p.save("evolution_coagulation_of_each_particle_kind.dat", clust_results)   



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


axes.plot(time,pos_results[:,6],"ro",\
          label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
# axes.plot(time,y_mean,"k+",\
#           label="<y>")
# axes.plot(time,z_mean,"mx",\
#           label="<z> ",linewidth=2.)
#p.ylim(1.4,2.6)
axes.legend()
p.xlabel('Time')
p.ylabel('delta position squared')
p.savefig("delta_r_squared.pdf")
 
p.clf()    






#print "pos_results_clusters_together is, ", pos_results_clusters_together


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


# axes.plot(time,pos_results_clusters_together[:,6],"ro",\
#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
# # axes.plot(time,y_mean,"k+",\
# #           label="<y>")
# # axes.plot(time,z_mean,"mx",\
# #           label="<z> ",linewidth=2.)
# #p.ylim(1.4,2.6)
# axes.legend()
# p.xlabel('Time')
# p.ylabel('delta position squared')
# p.savefig("delta_r_squared_clusters_together.pdf")
 
# p.clf()    






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


# axes.plot(time_selection,pos_results_selection[:,6],"ro",\
#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
# # axes.plot(time,y_mean,"k+",\
# #           label="<y>")
# # axes.plot(time,z_mean,"mx",\
# #           label="<z> ",linewidth=2.)
# #p.ylim(1.4,2.6)
# axes.legend()
# p.xlabel('Time')
# p.ylabel('delta position squared (selection)')
# p.savefig("delta_r_squared_selection.pdf")
 
# p.clf()    






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

axes.plot(time,pos_results[:,0],"ro",\
          label="<x>",linewidth=2.)
axes.plot(time,pos_results[:,1],"k+",\
          label="<y>")
axes.plot(time,pos_results[:,2],"mx",\
          label="<z> ",linewidth=2.)
#p.ylim(1.4,2.6)
#axes.legend()
p.xlabel('Time')
p.ylabel('Mean position')
p.savefig("mean_x_y_z.pdf")
 
p.clf()    






   

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


axes.plot(time,vel_results[:,6],"ro",\
          label="<v_x^2+v_y^2+v_z^2>",linewidth=2.)
# axes.plot(time,y_mean,"k+",\
#           label="<y>")
# axes.plot(time,z_mean,"mx",\
#           label="<z> ",linewidth=2.)
#p.ylim(1.4,2.6)
axes.legend()
p.xlabel('Time')
p.ylabel('delta velocity squared')
p.savefig("delta_v_squared.pdf")
 
p.clf()    





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

axes.plot(time,vel_results[:,0],"ro",\
          label="<v_x>",linewidth=2.)
axes.plot(time,vel_results[:,1],"k+",\
          label="<v_y>")
axes.plot(time,vel_results[:,2],"mx",\
          label="<v_z> ",linewidth=2.)
#p.ylim(1.4,2.6)
#axes.legend()
p.xlabel('Time')
p.ylabel('Mean velocity')
p.savefig("mean_vel_x_vel_y_vel_z.pdf")
 
p.clf()    



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


n_clus=s.sum(clust_results,axis=1)

print "the number of clusters is, ", n_clus


axes.plot(time,n_clus,"ro",\
          label="evolution total number of clusters")
types_arr=s.zeros(len(time))
types_arr[:]=types

axes.plot(time,types_arr,"b",\
          label="coagulation of all the individual clusters",linewidth=2.)


# axes.plot(time,y_mean,"k+",\
#           label="<y>")
# axes.plot(time,z_mean,"mx",\
#           label="<z> ",linewidth=2.)
#p.ylim(1.4,2.6)
axes.legend()
p.xlabel('Time')
p.ylabel('Number of clusters')
p.savefig("Coagulation_individual_clusters.pdf")
 
p.clf()    

#Now I think about a linear fit of delta_r_sq vs time

#         lin_fit=stats.linregress(s.log10(n_time_small),s.log10(r_time_small))
#         my_fra_low[count_config]=1./lin_fit[0]
#         r_stat_low[count_config]=lin_fit[2]

#I select the results only for the case when all the clusters are formed
#print "n_clus is, ", n_clus
#print "the shape of n_clus is, ", s.shape(n_clus)
cluster_formed=s.where(n_clus[:]==types)
#print "cluster_formed is, ", cluster_formed
#print "the shape of cluster_formed is, ", s.shape(cluster_formed)

time_formed=time[cluster_formed]
#print "time_formed is, ", time_formed
my_delta_r_sq=pos_results[:,6]
delta_r_sq_formed=my_delta_r_sq[cluster_formed]

p.save("delta_r_sq_time.dat",delta_r_sq_formed)

lin_fit=stats.linregress(time_formed,delta_r_sq_formed)

delta_r_sq_fitted=lin_fit[0]*time_formed+lin_fit[1]



Diff_coeff=lin_fit[0]

print "the linear coefficient [after cluster formation] is, ", lin_fit[0]




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


# axes.plot(time,pos_results[:,6],"ro",\
#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
# # axes.plot(time,y_mean,"k+",\
# #           label="<y>")
# # axes.plot(time,z_mean,"mx",\
# #           label="<z> ",linewidth=2.)
# #p.ylim(1.4,2.6)

# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)

# axes.legend()
# p.xlabel('Time')
# p.ylabel('delta position squared and fit')
# p.savefig("delta_r_squared_and_fit.pdf")
 
# p.clf()




# my_delta_r_sq=pos_results_clusters_together[:,6]
# delta_r_sq_formed=my_delta_r_sq[cluster_formed]


# lin_fit=stats.linregress(time_formed,delta_r_sq_formed)

# delta_r_sq_fitted=lin_fit[0]*time_formed+lin_fit[1]



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


# axes.plot(time,pos_results_clusters_together[:,6],"ro",\
#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
# # axes.plot(time,y_mean,"k+",\
# #           label="<y>")
# # axes.plot(time,z_mean,"mx",\
# #           label="<z> ",linewidth=2.)
# #p.ylim(1.4,2.6)

# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)

# axes.legend()
# p.xlabel('Time')
# p.ylabel('delta position squared and fit')
# p.savefig("delta_r_squared_and_fit_clusters_together.pdf")
 
# p.clf()    

# print "the linear coefficient is [clusters_together], ", lin_fit[0]


#Now the dynamics before the formation of the clusters


cluster_unformed=s.where(n_clus[:]!=types)
#print "cluster_formed is, ", cluster_formed
#print "the shape of cluster_formed is, ", s.shape(cluster_formed)

time_unformed=time[cluster_unformed]
#print "time_formed is, ", time_formed
#my_delta_r_sq=pos_results[:,6]
delta_r_sq_unformed=my_delta_r_sq[cluster_unformed]


lin_fit=stats.linregress(time_unformed,delta_r_sq_unformed)

delta_r_sq_fitted_unformed=lin_fit[0]*time_unformed+lin_fit[1]


print "the linear coefficient before forming the clusters is, ", lin_fit[0]



coag_choice=s.where((n_clus/n_clus[0])>0.2)



time_coag=time[coag_choice]
#print "time_formed is, ", time_formed
#my_delta_r_sq=pos_results[:,6]
delta_r_sq_coag=my_delta_r_sq[coag_choice]


print "delta_r_sq_coag is, ", delta_r_sq_coag


lin_fit=stats.linregress(time_coag,delta_r_sq_coag)

delta_r_sq_coag_fitted=lin_fit[0]*time_coag+lin_fit[1]




print "the linear coefficient at the early stages of coagulation is, ", lin_fit[0]








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


axes.plot(time,pos_results[:,6],"ro",\
          label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
# axes.plot(time,y_mean,"k+",\
#           label="<y>")
# axes.plot(time,z_mean,"mx",\
#           label="<z> ",linewidth=2.)
#p.ylim(1.4,2.6)

axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)

axes.plot(time_unformed, delta_r_sq_fitted_unformed,"k",linewidth=2.)

axes.plot(time_coag, delta_r_sq_coag_fitted,"m",linewidth=2.)



axes.legend()
p.xlabel('Time')
p.ylabel('delta position squared and fit')
p.savefig("delta_r_squared_and_fit.pdf")
 
p.clf()    




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


axes.plot(time,R_h_arr,"ro",\
          label="hydrodynamic radius",linewidth=2.)
# axes.plot(time,y_mean,"k+",\
#           label="<y>")
# axes.plot(time,z_mean,"mx",\
#           label="<z> ",linewidth=2.)
#p.ylim(1.4,2.6)

# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)

# axes.plot(time_unformed, delta_r_sq_fitted_unformed,"k",linewidth=2.)

# axes.plot(time_coag, delta_r_sq_coag_fitted,"m",linewidth=2.)



axes.legend()
p.xlabel('Time')
p.ylabel('R_h')
p.savefig("hydrodynamic_radius.pdf")
 
p.clf()    





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


axes.plot(time,R_gyr_arr,"ro",\
          label="gyration radius",linewidth=2.)
# axes.plot(time,y_mean,"k+",\
#           label="<y>")
# axes.plot(time,z_mean,"mx",\
#           label="<z> ",linewidth=2.)
#p.ylim(1.4,2.6)

# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)

# axes.plot(time_unformed, delta_r_sq_fitted_unformed,"k",linewidth=2.)

# axes.plot(time_coag, delta_r_sq_coag_fitted,"m",linewidth=2.)



axes.legend()
p.xlabel('Time')
p.ylabel('R_gyr')
p.savefig("gyration_radius.pdf")
 
p.clf()    




#print "R_h_arr is, ", R_h_arr

#Now I evaluate the ratio between the hydrodynamic radius and the radius of gyration as time goes by and I
# evaluate the 


R_h_over_R_gyr=s.mean(R_h_arr[cluster_formed]/R_gyr_arr[cluster_formed])


R_h_mean=R_h_arr[cluster_formed].mean()


R_gyr_mean=R_gyr_arr[cluster_formed].mean()


print "the ratio hydrodynamic radius over radius of gyration is, " ,R_h_over_R_gyr






#Now I take some time averages for a single cluster.

if (types==1): #i.e. I have a single cluster
    delta_r_sq_time_aver=s.var(pos_results[:,0:3], axis=0)
    delta_r_sq_time_aver=delta_r_sq_time_aver.sum()
    v_mean_time=s.mean(vel_results[:,0:3],axis=0)

    print "delta_r_sq_time_aver is, ", delta_r_sq_time_aver
    print "v_mean_time is, ", v_mean_time
    print "vel results are, ", vel_results

    print "the shape of vel_results is, ", s.shape(vel_results)


    print "<delta x squared>, ", pos_results[:,0].var()
    print "<delta y squared>, ", pos_results[:,1].var()
    print "<delta z squared>, ", pos_results[:,2].var()

    print "<x>, ", pos_results[:,0].mean()
    print "<y>, ", pos_results[:,1].mean()
    print "<z>, ", pos_results[:,2].mean()




    print "<v_x>, ", vel_results[:,0].mean()
    print "<v_y>, ", vel_results[:,1].mean()
    print "<v_z>, ", vel_results[:,2].mean()
    

    delta_v_sq_time_aver=s.var(vel_results[:,0:3], axis=0)
    delta_v_sq_time_aver=delta_v_sq_time_aver.sum()

    print "<delta v_x squared>, ", vel_results[:,0].var()
    print "<delta v_y squared>, ", vel_results[:,1].var()
    print "<delta v_z squared>, ", vel_results[:,2].var()
    

    

    print "delta_v_sq_time_aver is, ", delta_v_sq_time_aver









def Diff(kT,mu,r_eff):
    D=kT/(6.*s.pi*mu*r_eff)
    return D

def R_mob(kT,mu, D):
    r_eff=kT/(6.*s.pi*mu*D)

    return r_eff




kT=0.5
mu=1./(3.*s.pi)

Diff_R_h=Diff(kT,mu,R_h_mean)


Diff_R_gyr=Diff(kT,mu,R_gyr_mean)

print "D calculated with R_gyr is, ", Diff_R_gyr

print "D calculated with R_h is, ", Diff_R_h

R_mobility=R_mob(kT,mu,Diff_coeff)

print "The effective mobility radius is ", R_mobility

print "the ratio R_mobility/R_gyration is, ", R_mobility/R_gyr_mean

print "the ratio R_mobility/R_hydro_mean is, ", R_mobility/R_h_mean

# now I re-calculate the displacement for a given configuration as 









#Now I calculate the diffusion coefficient from the collected cluster positions

print "s.shape(R_cm_arr) is, ", s.shape(R_cm_arr)

#displacement=(R_cm_arr[:,0]-R_cm_arr)**2.


my_len=s.shape(R_cm_arr)[1]/3

p.save("centre_of_mass_positions.dat", R_cm_arr)


displacement=s.zeros((s.shape(R_cm_arr)[0],my_len))

print "R_cm_arr is, ", R_cm_arr[:,0:10]

for i in xrange(my_len):
    displacement[:,i]=s.sum(((R_cm_arr[:,0:3]-R_cm_arr[:,(i*3):((i+1)*3)])**2.),axis=1)

print "displacement is, ", displacement[:, 0:3]

p.save("displacement.dat", displacement)

# x_pos=x.zeros(my_len)
# y_pos=x.zeros(my_len)
# z_pos=x.zeros(my_len)

# for i in xrange(my_len):
#     x_pos=
    

displacement=displacement/6./time.max()

D_distr=displacement.mean(axis=1)

print "the trajectory-averaged distribution of diffusion coeff is, ", D_distr

D_mean=D_distr.mean()

print "the trajectory-averaged diffusion coeff is, ", D_mean

print "and its std is, ", D_distr.std()




#delta_r_sq_formed

i=70

print "I select time configuration (numbered from 0!!!), ", i

test_x_cm=R_cm_arr[:,(i*3)]
test_y_cm=R_cm_arr[:,(i*3+1)]
test_z_cm=R_cm_arr[:,(i*3+2)]
    
delta_x_sq_test=test_x_cm.var()
delta_y_sq_test=test_y_cm.var()
delta_z_sq_test=test_z_cm.var()

delta_r_sq_test=delta_x_sq_test+delta_y_sq_test+delta_z_sq_test

print "delta_r_sq_test is, ", delta_r_sq_test







print "So far so good"