← Back to branch summary

~jrjohansson/qutip/master

« back to all changes in this revision

Viewing changes to trilinear Monte-Carlo.py

  • Committer: Paul Nation
  • Date: 2011-04-21 04:46:56 UTC
  • Revision ID: git-v1:dd4c966b490aa468dfbd28cef66694df4bf235c8

Show diffs side-by-side

added added

removed removed

Lines of Context:
trilinear Monte-Carlo.py
 
1
from qutip import *
 
2
from pylab import *
 
3
import time
 
4
 
 
5
 
 
6
N0=15
 
7
N1=15
 
8
N2=15
 
9
K=1.0
 
10
gamma0=0.0
 
11
gamma1=0.0
 
12
gamma2=0.5
 
13
alpha=sqrt(5)
 
14
epsilon=0.5j #sqeezing parameter
 
15
tfinal=4.0
 
16
dt=0.05
 
17
tlist=arange(0.0,tfinal+dt,dt)
 
18
taulist=K*tlist #non-dimensional times
 
19
ntraj=50
 
20
 
 
21
#define operators
 
22
a0=tensor(destroy(N0),qeye(N1),qeye(N2))
 
23
a1=tensor(qeye(N0),destroy(N1),qeye(N2))
 
24
a2=tensor(qeye(N0),qeye(N1),destroy(N2))
 
25
 
 
26
num0=a0.dag()*a0
 
27
num1=a1.dag()*a1
 
28
num2=a2.dag()*a2
 
29
 
 
30
#dissipative operators for zero-temp. baths
 
31
C0=sqrt(2.0*gamma0)*a0
 
32
C1=sqrt(2.0*gamma1)*a1
 
33
C2=sqrt(2.0*gamma2)*a2
 
34
 
 
35
 
 
36
#initial state: coherent mode 0 & vacuum for modes #1 & #2
 
37
vacuum=tensor(basis(N0,0),basis(N1,0),basis(N2,0))
 
38
D=(alpha*a0.dag()-conj(alpha)*a0).expm()
 
39
psi0=D*vacuum
 
40
 
 
41
H=1j*K*(a0*a1.dag()*a2.dag()-a0.dag()*a1*a2)
 
42
Heff=H-0.5*1j*(C0.dag()*C0+C1.dag()*C1+C2.dag()*C2)
 
43
 
 
44
 
 
45
start_time=time.time()
 
46
ops=mcsolve(Heff,psi0,taulist,ntraj,[C0,C1,C2],[num0,num1,num2])
 
47
finish_time=time.time()
 
48
print 'time elapsed = ',finish_time-start_time
 
49
 
 
50
avg=sum(ops,axis=0)/ntraj
 
51
 
 
52
plot(taulist,avg[0],taulist,avg[1],taulist,avg[2])
 
53
show()
 
54