IQC/plotcurrent_mixed.py at master · ingstra/IQC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
from __future__ import division

import matplotlib.pylab as plt
import matplotlib.mlab as mlab
import numpy as np
from scipy.stats import norm
from scipy.misc import factorial
import sympy.mpmath as mp

current = np.loadtxt('current.dat')
h =  np.loadtxt('hermite.dat')
current=current[0,:]
dt = 1e-3
nruns = 10000

gamma = 2

mu = 0
variance = dt*nruns
sigma = np.sqrt(variance)

values, bins, _ = plt.hist(current,50,normed=1,histtype='bar')
area = sum(np.diff(bins)*values)
print area

#mufit, stdfit = norm.fit(current)
#print 'fitted mean: ',mufit #, 'calculated mean: '
#print 'fitted std: ',stdfit #,'calculated std: ',sigma

#plt.plot(h[:,0],h[:,1]**2,'.')

range = 5

x= np.linspace(-range,range,100)

n=1
h_poly = np.frompyfunc(mp.hermite,2,1)

const1 = 1/(np.sqrt(2**n * factorial(n)) * (2*np.pi)**0.25 )
rest1 = const1* np.exp(-x**2/4)
tot1  = (rest1*h_poly(n,x/np.sqrt(2)))**2

n=0
const0 = 1/(np.sqrt(2**n * factorial(n)) * (2*np.pi)**0.25 )
rest0 = const0* np.exp(-x**2/4)
tot0  = (rest0*h_poly(n,x/np.sqrt(2)))**2


plt.plot(x,tot1,'r-')


#plt.plot(traj[:,0],traj[:,1])
#plt.plot(exact[:,0],exact[:,1],'r--')

#plt.plot(traj[:,0],np.exp(-gamma*traj[:,0]),'r--')

#x = mu + sigma*np.random.randn(1000)
#n, bins, patches = plt.hist(x, 50, facecolor='green', alpha=0.75)
#y = mlab.normpdf( bins, mu, sigma)
#l = plt.plot(bins, y, 'r--')


#plt.plot(trace[:,0],trace[:,1])


def avg_Z(t):
    Omega = 1
    gamma = 1

    if 4*Omega < gamma:
        kappa = np.sqrt(np.power(gamma/4,2)-np.power(Omega,2))
    else:
        kappa = 1j*np.sqrt(-np.power(gamma/4,2)+np.power(Omega,2))

    K1 = np.power(Omega,2) / np.add(np.power(gamma,2) , 2*np.power(Omega,2) )
    hyp = np.cosh(kappa*t) + 3*gamma*np.sinh(kappa*t)/(4*kappa)

    return  K1*(1 - np.exp(-3*gamma*t/4)*hyp ) #- 1/2

def avg_Z_big_Omega(t):
    Omega = 5
    gamma = 1
    return -np.exp(-3*gamma*t/4)*np.cos(Omega*t)/2

def p_plus(t):
    Omega = 5
    gamma = 1
    return (1 - np.exp(-3*gamma*t/4)*np.cos(Omega*t))/2

#plt.plot(traj[:,0],avg_Z(traj[:,0]),'r--')

#plt.plot(data[:,0],p_plus(data[:,0]),'g')


plt.show()