← Back to branch summary

~jrjohansson/qutip/master

« back to all changes in this revision

Viewing changes to test_ode2es.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:
test_ode2es.py
 
1
from qutip import *
 
2
from pylab import *
 
3
 
 
4
def probevolve(E,kappa,gamma,g,wc,w0,wl,N,tlist):
 
5
 
 
6
    ida    = qeye(N)
 
7
    idatom = qeye(2)
 
8
 
 
9
    # Define cavity field and atomic operators
 
10
    a  = tensor(destroy(N),idatom)
 
11
    sm = tensor(ida,sigmam())
 
12
 
 
13
    # Hamiltonian
 
14
    H = (w0-wl)*sm.dag()*sm + (wc-wl)*a.dag()*a + 1j*g*(a.dag()*sm - sm.dag()*a) + E*(a.dag()+a)
 
15
 
 
16
    #collapse operators
 
17
    C1=sqrt(2*kappa)*a
 
18
    C2=sqrt(gamma)*sm
 
19
    C1dC1=C1.dag()*C1
 
20
    C2dC2=C2.dag()*C2
 
21
 
 
22
    #intial state
 
23
    psi0 = tensor(basis(N,0),basis(2,1))
 
24
    rho0 = psi0 * trans(psi0);
 
25
 
 
26
    # Calculate the Liouvillian
 
27
    c_op_list = [C1, C2]
 
28
    n_op      = len(c_op_list)
 
29
    L = -1j * (spre(H) - spost(H))
 
30
    for m in range(0, n_op):
 
31
        cdc = c_op_list[m].dag() * c_op_list[m]
 
32
        L += spre(c_op_list[m])*spost(c_op_list[m].dag())-0.5*spre(cdc)-0.5*spost(cdc)
 
33
 
 
34
       
 
35
    # Calculate solution as an exponential series
 
36
    start_time=time.time()
 
37
    rhoES = ode2es(L,rho0);
 
38
    print 'time elapsed (ode2es) = ' +str(time.time()-start_time) 
 
39
    
 
40
    # Calculate expectation values
 
41
    start_time=time.time()
 
42
    count1  = esval(scalar_expect(C1dC1,rhoES),tlist);
 
43
    count2  = esval(scalar_expect(C2dC2,rhoES),tlist);
 
44
    infield = esval(scalar_expect(a,rhoES),tlist);
 
45
    print 'time elapsed (esval) = ' +str(time.time()-start_time) 
 
46
 
 
47
    return count1, count2, infield
 
48
 
 
49
 
 
50
#-----------------------------------------------------------------------------
 
51
# 
 
52
#--------------------------------------------------------------------------
 
53
kappa = 2; 
 
54
gamma = 0.2;
 
55
g  = 1; 
 
56
wc = 0; 
 
57
w0 = 0; 
 
58
wl = 0;
 
59
N  = 4; 
 
60
E  = 0.5;
 
61
 
 
62
tlist = linspace(0,10,200);
 
63
 
 
64
start_time=time.time()
 
65
[count1,count2,infield] = probevolve(E,kappa,gamma,g,wc,w0,wl,N,tlist);
 
66
print 'time elapsed = ' +str(time.time()-start_time) 
 
67
 
 
68
plot(tlist,real(count1))
 
69
plot(tlist,real(count2))
 
70
xlabel('Time')
 
71
ylabel('Transmitted Intensity and Spontaneous Emission')
 
72
show()
 
73