← Back to branch summary

~jrjohansson/qutip/master

« back to all changes in this revision

Viewing changes to test_corr4.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_corr4.py
 
1
 
 
2
from qutip import *
 
3
from pylab import *
 
4
 
 
5
#
 
6
#
 
7
#
 
8
def probcorr(E,kappa,gamma,g,wc,w0,wl,N):
 
9
    #
 
10
    # [corrES,covES] = probcorr(E,kappa,gamma,g,wc,w0,wl,N)
 
11
    #  returns the two-time correlation and covariance of the intracavity 
 
12
    #  field as exponential series for the problem of a coherently driven 
 
13
    #  cavity with a two-level atom
 
14
    #
 
15
    #  E = amplitude of driving field, kappa = mirror coupling,
 
16
    #  gamma = spontaneous emission rate, g = atom-field coupling,
 
17
    #  wc = cavity frequency, w0 = atomic frequency, wl = driving field frequency,
 
18
    #  N = size of Hilbert space for intracavity field (zero to N-1 photons)
 
19
    #  
 
20
 
 
21
    ida    = qeye(N)
 
22
    idatom = qeye(2)
 
23
 
 
24
    # Define cavity field and atomic operators
 
25
    a  = tensor(destroy(N),idatom)
 
26
    sm = tensor(ida,sigmam())
 
27
 
 
28
    # Hamiltonian
 
29
    H = (w0-wl)*sm.dag()*sm + (wc-wl)*a.dag()*a + 1j*g*(a.dag()*sm - sm.dag()*a) + E*(a.dag()+a)
 
30
 
 
31
    #collapse operators
 
32
    C1=sqrt(2*kappa)*a
 
33
    C2=sqrt(gamma)*sm
 
34
 
 
35
    # Calculate the Liouvillian
 
36
    c_op_list = [C1, C2]
 
37
    n_op      = len(c_op_list)
 
38
    L = -1.0j * (spre(H) - spost(H))
 
39
    for m in range(0, n_op):
 
40
        cdc = c_op_list[m].dag() * c_op_list[m]
 
41
        L += spre(c_op_list[m])*spost(c_op_list[m].dag())-0.5*spre(cdc)-0.5*spost(cdc)
 
42
 
 
43
    # Find steady state density matrix and field
 
44
    rhoss = steady(L);
 
45
   
 
46
    # Initial condition for regression theorem
 
47
    arhoad = a * rhoss * a.dag();
 
48
    #psi = tensor(basis(N,0),basis(2,0))
 
49
    #arhoad = psi * trans(psi)
 
50
    #arhoad = rhoss
 
51
    #print "rhoss =", rhoss
 
52
 
 
53
    # Solve differential equation with this initial condition
 
54
    solES = ode2es(L, arhoad);
 
55
 
 
56
    # Find trace(a' * a * solution)
 
57
    corrES = scalar_expect((a.dag() * a), solES);
 
58
 
 
59
    return corrES
 
60
 
 
61
#
 
62
#
 
63
#
 
64
kappa = 2; 
 
65
gamma = 0.2; 
 
66
g = 5;
 
67
wc = 0; 
 
68
w0 = 0; 
 
69
wl = 0;
 
70
N = 5; 
 
71
E = 0.5;
 
72
 
 
73
start_time=time.time()
 
74
corrES = probcorr(E,kappa,gamma,g,wc,w0,wl,N);
 
75
print 'time elapsed (probcorr) = ' +str(time.time()-start_time) 
 
76
 
 
77
tlist = linspace(0,10.0,200);
 
78
 
 
79
start_time=time.time()
 
80
corr  = esval(corrES,tlist);
 
81
print 'time elapsed (esval) = ' +str(time.time()-start_time) 
 
82
 
 
83
figure(1)
 
84
plot(tlist,real(corr))
 
85
xlabel('Time')
 
86
ylabel('Correlation and covariance')
 
87
show()
 
88
 
 
89