myprojects-hg-reborn

Fork

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

from scipy.special import * # to be able to use the error function

def vth(q,beta):
	vth=sqrt(q/(2./beta))
	return vth

def discr(alpha,beta):
	discr=sqrt(beta**2.+4.*alpha*beta)
	return discr


def lambda1(alpha,beta): # do not reintroduce alpha and beta! they are already in discriminant!
	lambda1= -0.5*(beta+discr(alpha,beta))
	return lambda1
	


def lambda2(alpha,beta):
	lambda2=0.5*(-beta+discr(alpha,beta))
	return lambda2




def aux0(alpha,beta):
	aux0=lambda1(alpha,beta)-lambda2(alpha,beta)
	return aux0



def aux2(alpha,beta):
	aux2=lambda1(alpha,beta)*lambda2(alpha,beta)
	return aux2

def aux3(alpha,beta):
	aux3=lambda1(alpha,beta)+lambda2(alpha,beta)
	return aux3



#Time propagators

def GXX(alpha,beta,dT):
	GXX=(lambda1(alpha,beta)*exp(lambda2(alpha,beta)*dT)-lambda2(alpha,beta) \
	*exp(lambda1(alpha,beta)*dT)) \
	/aux0(alpha,beta)
	return GXX


def GXVx(alpha,beta,dT):
	GXVx=(exp(lambda1(alpha,beta)*dT)-exp(lambda2(alpha,beta)*dT))/aux0(alpha,beta)
	return GXVx


def GVxX(alpha,beta,dT):
	GVxX=aux2(alpha,beta)*(exp(lambda2(alpha,beta)*dT)-exp(lambda1(alpha,beta)*dT)) \
	/aux0(alpha,beta)
	return GVxX

def GVxVx(alpha,beta,dT):
	GVxVx=(lambda1(alpha,beta)*exp(lambda1(alpha,beta)*dT)- \
	lambda2(alpha,beta)*exp(lambda2(alpha,beta)*dT))/aux0(alpha,beta)
	return GVxVx

#variances

def varX(q,alpha,beta,dT):
	if alpha !=0.:
		varX=q*(0.5*(exp(2.*lambda1(alpha,beta)*dT)-1.)/lambda1(alpha,beta)-  \
		2.*(exp(aux3(alpha,beta)*dT)-1.)/aux3(alpha,beta) \
		+0.5*(exp(2.*lambda2(alpha,beta)*dT)-1.)/lambda2(alpha,beta)) \
		/aux0(alpha,beta)**2.
		return varX
	else:
		varX=q*(2*beta*dT-3.+4.*exp(-beta*dT)-exp(-2.*beta*dT))/(2.*beta**3.)
		return varX
	
	

def varVx(q,alpha,beta,dT):
	varVx=q*(0.5*(exp(2.*lambda1(alpha,beta)*dT)-1.)*lambda1(alpha,beta)- \
	2.*aux2(alpha,beta)*(exp(aux3(alpha,beta))-1.)/aux3(alpha,beta)+ \
	0.5*(exp(2.*lambda2(alpha,beta)*dT)-1.)*lambda2(alpha,beta))/aux0(alpha,beta)**2.
	return varVx

def covXVx(q,alpha,beta,dT):
	covXVx=q*0.5*(exp(lambda1(alpha,beta)*dT)-exp(lambda2(alpha,beta)*dT))**2. \
	/aux0(alpha,beta)**2.
	return covXVx

vec_covxvx=vectorize(covXVx)

#mean values

def meanX(q,alpha,beta,dT,x0,vx0):
	meanx=GXX(alpha,beta,dT)*x0+GXVx(alpha,beta,dT)*vx0
	return meanx

vec_meanx=vectorize(meanX)

def meanVx(q,alpha,beta,dT,x0,vx0):
	meanVx=GVxX(alpha,beta,dT)*x0+GVxVx(alpha,beta,dT)*vx0
	return meanVx


vec_meanvx=vectorize(meanVx)

# now the determinant

def detSigma(q,alpha,beta,dT):
	deter=varX(q,alpha,beta,dT)*varVx(q,alpha,beta,dT) \
	-(covXVx(q,alpha,beta,dT))**2.
	return deter

vec_detsigma=vectorize(detSigma)



# now the part dealing with the propagator

def tinv(q,alpha,beta,dT):
	tinv=varX(q,alpha,beta,dT)/(2.*detSigma(q,alpha,beta,dT))
	return tinv

def Veff(q,alpha,beta,dT,x,vx):
	Veff=vx-covXVx(q,alpha,beta,dT)*x/varX(q,alpha,beta,dT)
	return Veff

def V0eff(q,alpha,beta,dT,x0,vx0):
	V0eff=meanVx(q,alpha,beta,dT,x0,vx0)-covXVx(q,alpha,beta,dT)* \
	meanX(q,alpha,beta,dT,x0,vx0)/varX(q,alpha,beta,dT)
	return V0eff

def Xeff(q,alpha,beta,dT,x):
	Xeff=sqrt(detSigma(q,alpha,beta,dT))*x/abs(varX(q,alpha,beta,dT))
	return Xeff

def X0eff(q,alpha,beta,dT,x0,vx0):
	X0eff=sqrt(detSigma(q,alpha,beta,dT))*meanX(q,alpha,beta,dT,x0,vx0) \
	/abs(varX(q,alpha,beta,dT))
	return X0eff

def r(q,alpha,beta,dT,x,vx):
	r=sqrt(Veff(q,alpha,beta,dT,x,vx)**2.+Xeff(q,alpha,beta,dT,x)**2.)
	return r


def rprime(q,alpha,beta,dT,x0,vx0):
	rprime=sqrt(V0eff(q,alpha,beta,dT,x0,vx0)**2.+X0eff(q,alpha,beta,dT,x0,vx0)**2.)
	return rprime



def theta(q,alpha,beta,dT,x,vx):
	theta=(1.-sign(Xeff(q,alpha,beta,dT,x)/Veff(q,alpha,beta,dT,x,vx)))*pi/2.+ \
	arctan(Xeff(q,alpha,beta,dT,x)/Veff(q,alpha,beta,dT,x,vx))+pi
	return theta


def thetaprime(q,alpha,beta,dT,x0,vx0):
	thetaprime=(1.-sign(X0eff(q,alpha,beta,dT,x0,vx0)/V0eff(q,alpha,beta,dT,x0,vx0)))\
	*pi/2+ arctan(X0eff(q,alpha,beta,dT,x0,vx0)/V0eff(q,alpha,beta,dT,x0,vx0))+pi
	return thetaprime

def a(q,alpha,beta,dT,x,vx,x0,vx0):
	a=2.*sqrt(r(q,alpha,beta,dT,x,vx)*rprime(q,alpha,beta,dT,x0,vx0)\
	*tinv(q,alpha,beta,dT))\
	*cos((theta(q,alpha,beta,dT,x,vx)-thetaprime(q,alpha,beta,dT,x0,vx0))/2.)
	return a



def R2(q,alpha,beta,dT,x,vx,x0,vx0):
	R2=r(q,alpha,beta,dT,x,vx)**2.-2.*rprime(q,alpha,beta,dT,x0,vx0)\
	*cos(theta(q,alpha,beta,dT,x,vx)-thetaprime(q,alpha,beta,dT,x0,vx0))
	return R2

def MultiGreenProp(q,alpha,beta,dT,x,vx,x0,vx0):
	Green=exp(-R2(q,alpha,beta,dT,x,vx,x0,vx0)*tinv(q,alpha,beta,dT))\
	*(1.-erf(-a(q,alpha,beta,dT,x,vx,x0,vx0)))\
	/(4.*pi*sqrt(detSigma(q,alpha,beta,dT)))
	return Green

# now the antipropagator
def thetaprimeAnti(q,alpha,beta,dT,x0,vx0):
	thetaAnti=-thetaprime(q,alpha,beta,dT,x0,vx0)
	return thetaAnti

def aAnti(q,alpha,beta,dT,x,vx,x0,vx0):
	aAnti=2.*sqrt(r(q,alpha,beta,dT,x,vx)*rprime(q,alpha,beta,dT,x0,vx0)\
	*tinv(q,alpha,beta,dT))\
	*cos((theta(q,alpha,beta,dT,x,vx)-thetaprimeAnti(q,alpha,beta,dT,x0,vx0))/2.)
	return aAnti

	
def R2Anti(q,alpha,beta,dT,x,vx,x0,vx0):
	R2Anti=r(q,alpha,beta,dT,x,vx)**2.+rprime(q,alpha,beta,dT,x0,vx0)**2.\
	-2.*r(q,alpha,beta,dT,x,vx)*rprime(q,alpha,beta,dT,x0,vx0)\
	*cos(theta(q,alpha,beta,dT,x,vx)-thetaprimeAnti(q,alpha,beta,dT,x0,vx0))
	return R2Anti

def MultiGreenAnti(q,alpha,beta,dT,x,vx,x0,vx0):
	GreenAnti=exp(-R2(q,alpha,beta,dT,x,vx,x0,vx0)*tinv(q,alpha,beta,dT))\
	*(1.-erf(-aAnti(q,alpha,beta,dT,x,vx,x0,vx0)))\
	/(4.*pi*sqrt(detSigma(q,alpha,beta,dT)))
	return GreenAnti


def Green_complete(q,alpha,beta,dT,x,vx,x0,vx0):
	G_compl=MultiGreenProp(q,alpha,beta,dT,x,vx,x0,vx0)\
	-MultiGreenAnti(q,alpha,beta,dT,x,vx,x0,vx0)
	return G_compl

vec_green=vectorize(Green_complete)

alpha=0.
beta=1.
q=2.
vx0=-1.
x0=20.
x=1e-10
vx=linspace(-8.,8.,2000)
T=9000.

my_green=vec_green(q,alpha,beta,T,x,vx,x0,vx0)

a1=aux0(alpha,beta)
a2=aux2(alpha,beta)
a3=aux3(alpha,beta)


z=aux0(alpha,beta)
print 'aux0 is', z
z=covXVx(q,alpha,beta,1.)
print 'covXvx is', z
z=detSigma(q,alpha,beta,9.)

print 'detsigma is', z

dT=linspace(1e-3,1e1,100)
my_meanx=vec_meanx(q,alpha,beta,dT,x0,vx0)
my_meanvx=vec_meanvx(q,alpha,beta,dT,x0,vx0)
my_detsigma=vec_detsigma(q,alpha,beta,dT)
covariance=vec_covxvx(q,alpha,beta,dT)

print 'my cov is', covariance

#print 'my_meanvx is', my_meanvx
#print 'a1,a2,a3 are',a1,a2,a3



pylab.plot(dT,my_meanx)
pylab.xlabel('Time')
pylab.ylabel('Expectation value of X')
pylab.title('Mean X')
pylab.grid(True)
pylab.savefig('mean_X')
#pylab.show() # uncomment them to plot a distribution

pylab.hold(False)

pylab.plot(dT,my_meanvx)
pylab.xlabel('Time')
#pylab.ylabel('Expectation value of Vx')
pylab.title('Mean Vx')
pylab.grid(True)
pylab.savefig('mean_Vx')


pylab.plot(dT,covariance)
pylab.xlabel('Time')
#pylab.ylabel('Expectation value of Vx')
pylab.title('covxvx')
pylab.grid(True)
pylab.savefig('covariance')


pylab.plot(dT,my_detsigma)
pylab.xlabel('Time')
#pylab.ylabel('Expectation value of Vx')
pylab.title('determinant')
pylab.grid(True)
pylab.savefig('determinant')

#pylab.hold(True)

my_result=pylab.load("wall_T500.txt")

my_result2=pylab.load("wall_T9000.txt")

def asym_v(v):
	if (v>=0):
		v_prof=sqrt(abs(v))*exp((-v**2.)/2.)
		return v_prof
	elif (v<0):
		v_prof=0.
		return v_prof

vec_v_prof=vectorize(asym_v)

v_asym=vec_v_prof(vx)

pylab.plot(vx,my_green/max(my_green),my_result[:,0],my_result[:,1]/max(my_result[:,1])\
,my_result2[:,0],my_result2[:,1]/max(my_result2[:,1]),linewidth=1.5)
pylab.legend(('analytics','T=500','T=9000'))
pylab.xlabel('velocity')
#pylab.ylabel('Expectation value of Vx')
pylab.title('Green function')
pylab.grid(True)
pylab.savefig('green_function_wall')
#pylab.show()

pylab.hold(False)

pylab.plot(my_result2[:,0],my_result2[:,1]/max(my_result2[:,1])\
,vx,v_asym/max(v_asym),linewidth=1.5)
pylab.legend(('T=9000', 'asymptotics'))
pylab.xlabel('velocity')
pylab.ylabel('Expectation value of Vx')
pylab.title('Green function')
pylab.grid(True)
pylab.savefig('green_function_wall-comparison')
#pylab.show()

def find_max(my_arr):
	for i in range(0,len(my_arr)):
		if my_arr[i]==max(my_arr):
			max_pos=i
			return max_pos

m_anal=find_max(v_asym)
print 'the maximum of v_asym is at',vx[m_anal]


m_num=find_max(my_result2[:,1])
print 'the maximum of the numerical result is at',my_result2[m_num,0]

print 'so far so good'