提交: 8b4222acf1f0c5903aa23a88c6e32bf4f834ced2 - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修订版: 8b4222acf1f0c5903aa23a88c6e32bf4f834ced2 (tree)
• 时间: 2007-12-06 23:41:32
• 作者: iselllo
• Commiter: iselllo

Log Message

I updated plot_statistics.py; now the code is also able to calculate the
velocity autocorrelation function.

更改概述

• modified: Python-codes/plot_statistics.py

差异

diff -r cd60708cc67b -r 8b4222acf1f0 Python-codes/plot_statistics.py
--- a/Python-codes/plot_statistics.py	Tue Dec 04 12:34:22 2007 +0000
+++ b/Python-codes/plot_statistics.py	Thu Dec 06 14:41:32 2007 +0000
@@ -5,11 +5,31 @@
 import pylab as p
 
 n_part=3000
-n_config=100
+n_config=100 # total number of configurations, which are numbered from 0 to n_config-1
 t_step=1.
 
 config_ini=50
 
+time=s.linspace(0.,((n_config-1)*t_step),n_config) #I define the time along the simulation
+
+
+
+
+#some physical parameters of the system
+T=0.1  #temperature
+beta=0.1 #damping coefficient
+v_0=0.  #initial particle mean velocity
+
+
+n_s_ini=40  #element in the time-sequence corresponding
+ # to the chosen reference time for calculating the velocity autocorrelation
+ #function
+
+s_ini=(n_s_ini)*t_step #I define a specific time I will need for the calculation
+# of the autocorrelation function
+
+print 'the time chosen for the velocity correlation function is, ', s_ini
+print 'and the corresponding time on the array is, ', time[n_s_ini]
 
 
 
@@ -60,14 +80,17 @@
 	mean_x_arr=s.zeros(n_config)
 	var_x_arr=s.zeros(n_config)
 
-	for i in xrange(0, n_config):
-		mean_x_arr[i]=s.mean(x_arr[i,:])
-		var_x_arr[i]=s.var(x_arr[i,:])
+	#for i in xrange(0, n_config):
+		#mean_x_arr[i]=s.mean(x_arr[i,:])
+		#var_x_arr[i]=s.var(x_arr[i,:])
+
+	mean_x_arr[:]=s.mean(x_arr,axis=1)
+	var_x_arr[:]=s.var(x_arr,axis=1)
+
 
 	print "mean_x_arr is, ", mean_x_arr
 	print "var_x_arr is, ", var_x_arr
 
-	time=s.linspace(0.,(n_config*t_step),n_config)
 
 	p.plot(time, var_x_arr ,linewidth=2.)
 	p.xlabel('time')
@@ -79,24 +102,105 @@
 	p.hold(False)
 
 	p.save("mean_square_disp.dat", var_x_arr)
+# Now I do pretty much the same thing with the velocity
 
+	tot_config_vel=p.load("total_configuration_vel.dat")
+
+	
+	 
+	  
+	   
+	print "the length of tot_config is, ", len(tot_config)
+	tot_config_vel=s.reshape(tot_config_vel,(n_config,3*n_part))
+
+	#print "tot_config at line 10 is, ", tot_config[10,:]
+
+
+	v_0=tot_config_vel[0,0] #initial x position of all my particles
+	v_arr=s.zeros((n_config,n_part))
+	print 'OK creating v_arr'
+	v_col=s.arange(n_part)*3
+
+	#print "x_col is, ", x_col
+
+	for i in xrange(0,n_part):
+		v_arr[:,i]=tot_config_vel[:,3*i]
+
+	mean_v_arr=s.zeros(n_config)
+	var_v_arr=s.zeros(n_config)
+
+#	for i in xrange(0, n_config):
+	#	mean_v_arr[i]=s.mean(v_arr[i,:])
+		#var_v_arr[i]=s.var(v_arr[i,:])
+
+	mean_v_arr[:]=s.mean(v_arr,axis=1)
+	var_v_arr[:]=s.var(v_arr,axis=1)
+
+
+
+
+	#print "mean_v_arr is, ", mean_v_arr
+	#print "var_v_arr is, ", var_v_arr
+
+	#time=s.linspace(0.,(n_config*t_step),n_config)
+
+	p.plot(time, var_v_arr ,linewidth=2.)
+	p.xlabel('time')
+	p.ylabel('velocity variance')
+	#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
+	p.title('VAR(v) vs time')
+	p.grid(True)
+	p.savefig('velocity_variance.pdf')
+	p.hold(False)
+
+	p.save("velocity_variance.dat", var_v_arr)
+
+	#Now I want to calculate the autocorrelation function.
+	#First, I need to fix an intial time. I choose something which is well past the
+	#ballistic regime
+	
+
+
+	#Now I start calculating the velocity autocorrelation function
+	
+	vel_autocor=s.zeros(n_config)
+
+	
+	
+       	for i in xrange(0, n_config):
+		vel_autocor[i]=s.mean(v_arr[i]*v_arr[n_s_ini])
+
+
+	p.plot(time, vel_autocor ,linewidth=2.)
+	p.xlabel('time')
+	p.ylabel('<v(s)v(t)> ')
+	#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
+	p.title('Velocity correlation')
+	p.grid(True)
+	p.savefig('velocity_correlation_numerics.pdf')
+	p.hold(False)
+
+
+	p.save("velocity_autocorrelation.dat", vel_autocor)
+
+
+		
+
+	
 elif (calculate ==0):
 	
 	var_x_arr=p.load("mean_square_disp.dat")
 	
 	#now some analytics of the mean-square displacement
-	def x_sq(T,Gamma,v_0 ,t):
-		z=2.*T/(Gamma)*t+v_0**2./Gamma**2.*(1.-s.exp(-Gamma*t))\
-		-T/Gamma**2.*(3.-4.*s.exp(-Gamma*t)+s.exp(-2.*Gamma*t))
+	def x_sq(T,beta,v_0 ,t):
+		z=2.*T/(beta)*t+v_0**2./beta**2.*(1.-s.exp(-beta*t))\
+		-T/beta**2.*(3.-4.*s.exp(-beta*t)+s.exp(-2.*beta*t))
 		return z
 	
-	time=s.linspace(0.,(n_config*t_step),n_config)
+#	time=s.linspace(0.,(n_config*t_step),n_config)
 	
-	T=0.1
-	Gamma=0.1
-	v_0=0.
 	
-	x_sq_book=x_sq(T,Gamma,v_0,time)
+	x_sq_book=x_sq(T,beta,v_0,time)
 	
 	p.plot(time, var_x_arr,time, x_sq_book ,linewidth=2.)
 	p.xlabel('time')
@@ -106,8 +210,84 @@
 	p.grid(True)
 	p.savefig('mean_square_disp_comparison_analytics.pdf')
 	p.hold(False)
+	
+	#now a similar thing for the velocity
+	
+	def v_sq(T,beta,v_0,Gamma,t):
+		z=Gamma**2./(2.*beta)+(v_0**2.-Gamma**2./(2.*beta))*s.exp(-2.*beta*t)
+		return z
+	
+	Gamma=s.sqrt(2.*beta*T)
+	
+	#NB: Gamma (NOT the damping parameter which I call beta here!!!) is related to the nosie strength i.e.
+	# Gamma**2./(2*beta)=k_BT/m i.e. Gamma=sqrt(2*beta*k_b*T/m), but I work in units where
+	#Gamma=sqrt(2*beta*T)
+	
+	
+	v_sq_book=v_sq(T,beta,v_0,Gamma,time)
+	
+	var_v_arr=p.load("velocity_variance.dat")
+	
+	
+	p.plot(time, var_v_arr,time, v_sq_book ,linewidth=2.)
+	p.xlabel('time')
+	p.ylabel('velocity variance ')
+	#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
+	p.title('VAR(x) vs time')
+	p.grid(True)
+	p.savefig('velocity_variance_comparison_analytics.pdf')
+	p.hold(False)
+
+	def vel_corr(t,s_ini,v_0,beta,T):
+		z=(v_0**2.-T)*s.exp(-beta*(t+s_ini))+T*s.exp(-beta*abs(t-s_ini))
+#		z=(v_0**2.-T)*s.exp(-t/10.)
+		return z
+
+#	def vel_corr_asym(t,s,beta,T):
+	#	z=T*s.exp(-beta*abs(t-s))
+		#return z
+
+
+
+
+	print 'the initial time for the velocity correlation is',s_ini
+
+	print 'T is,', T
+	print 'beta is, ', beta
+	print 'time is, ', time
+	print 'v_0 is, ', v_0
+	my_vel_corr= vel_corr(time,s_ini,v_0,beta,T)
+	print 'my_vel_corr is, ', my_vel_corr
+
 
 
 	
 
-print "So far so good"
\ No newline at end of file
+	print 'OK before plot'
+	
+	p.plot(time, my_vel_corr ,linewidth=2.)
+	p.xlabel('time')
+	p.ylabel('<v(s)v(t)> ')
+	#p.legend(('beta=1e-2,100 part','beta=1e-1, 100 part', 'beta=1e-1, 200 part'))
+	p.title('Velocity correlation')
+	p.grid(True)
+	p.savefig('velocity_correlation_analytics.pdf')
+	p.hold(False)
+
+	#Now I read the velocity autocorrelation function I calculated numerically
+	# and try seeing how it compares with the analytics
+
+	vel_autocor_numerics=p.load("velocity_autocorrelation.dat")
+	
+
+       	p.plot(time, vel_autocor_numerics,time, my_vel_corr ,linewidth=2.)
+	p.xlabel('time')
+	p.ylabel('velocity autocorrelation ')
+	p.legend(('Numerical Result','Analytical Result'))
+	p.title('Velocity Autocorrelation')
+	p.grid(True)
+	p.savefig('velocity_autocorrelation_comparison_numerics_analytics.pdf')
+	p.hold(False)
+
+
+print "So far so good"