initial addition

1 files changed, 190 insertions, 0 deletions

diff --git a/liouville.m b/liouville.m
new file mode 100644
index 0000000..baff45a
--- /dev/null
+++ b/liouville.m
@@ -0,0 +1,190 @@
+
+levels=1:3;
+
+%      -----------  |1>
+%         /     \
+%    E_d /       \ 
+%       /         \ E_p
+%      /           \
+%    -------- |3>   \
+%                    \
+%               ___________ |2>
+
+
+% some physical constants
+hbar=1;
+im_one=0+1i;
+
+Nlevels=3;
+w12=1e9;
+w_hpf=6800;
+w13=w12-w_hpf
+
+
+% unperturbed Hamiltonian energy levels
+levels_energy=[ w12, 0,  w_hpf];
+levels_energy=levels_energy*hbar;
+H0=zeros(Nlevels);
+H0=diag(levels_energy);
+%for i=1:Nlevels
+	%H0(i,i)=levels_energy(i);
+%endfor
+
+% decay matrix g(i,j) correspnds to decay from i-->j
+gamma=6;
+gamma_13=.01
+g_decay=zeros(Nlevels);
+g_decay(1,2)=gamma; %upper level decay
+g_decay(1,3)=gamma; %upper level decay
+g_decay(3,2)=gamma_13; % lower levels mixing
+g_decay(2,3)=gamma_13; % lower levels mixing
+
+%defasing matris
+g_deph=0;
+g_dephasing=zeros(Nlevels);
+g_dephasing(1,2)=g_deph;
+g_dephasing(2,1)=g_dephasing(1,2);
+g_dephasing(1,3)=g_deph;
+g_dephasing(3,1)=g_dephasing(1,3);
+
+
+
+% dipole matrix
+dipole_elements=zeros(Nlevels);
+dipole_elements(1,2)=1;
+dipole_elements(2,1)=dipole_elements(1,2);
+dipole_elements(3,1)=1;
+dipole_elements(1,3)=dipole_elements(3,1);
+
+
+%EM field definition
+Ep=100.0;
+Epc=conj(Ep);
+Ed=0;
+Edc=conj(Ed);
+wd=w13;
+wp=w12;
+modulation_freq=[0,  wp, wd, -wp, -wd, wp-wd, wd-wp];
+E_field        =[0,  Ep, Ed, Epc, Edc, 0,     0    ];
+Nfreq=length(modulation_freq);
+
+% now we create Liouville indexes list
+% we unwrap density matrix and assign all posible
+% frequencies as well
+% resulting vector should be Nlevels x Nlevels x length(modulation_freq)
+i=0;
+for w=1:length(modulation_freq)
+	for r=1:Nlevels
+		for c=1:Nlevels
+			i+=1;
+			rhoLiouville(i)=0;
+			rhoLiouville_w(i)=w;
+			rhoLiouville_r(i)=r;
+			rhoLiouville_c(i)=c;
+		endfor
+	endfor
+endfor
+
+\
+
+N=length(rhoLiouville);
+% Liouville operator matrix 
+L=zeros(N); % NxN matrix
+Li=zeros(N); % NxN Liouville interactive
+L0=zeros(N); % NxN Liouville from unperturbed hamiltonian
+
+function ret=decay_total(g_decay,i)
+	ret=0;
+	for k=1:size(g_decay)(1)
+		ret=ret+g_decay(i,k);
+	endfor
+endfunction
+
+function ret=kron_delta(i,j)
+	ret=((i==j));
+endfunction
+
+for p=1:N
+	for s=1:N
+		j=rhoLiouville_r(p);
+		k=rhoLiouville_c(p);
+		m=rhoLiouville_r(s);
+		n=rhoLiouville_c(s);
+		% we garanted to know frequency of final and initial rhoLiouville
+		w1i=rhoLiouville_w(p);
+		w2i=rhoLiouville_w(s);
+		w_jk=modulation_freq(w1i);
+		w_mn=modulation_freq(w2i);
+		% thus we know L matrix element frequency which we need to match
+		w_l=w_jk-w_mn;
+		% lets search this wrequency in the list of available frequencyes
+		% but since we not garanteed to find it lets assign temporary 0 to Liouville matrix element
+		L(p,s)=0;
+		for w3i=1:length(modulation_freq)
+			w_iner=modulation_freq(w3i);
+			if ((w_iner == w_l))
+				%such frequency exist in the list of modulation frequencies
+				if ((w_iner == 0))
+					L0=H0(j,m)*kron_delta(k,n)-H0(n,k)*kron_delta(j,m);
+					decay_part=\
+						( decay_total(g_decay,k)/2 + decay_total(g_decay,j)/2 + g_dephasing(j,k) )* kron_delta(j,m)*kron_delta(k,n) \
+						- kron_delta(m,n)*kron_delta(j,k)*g_decay(m,j) ;
+					decay_part=decay_part*hbar/im_one;
+					L0=decay_part;
+					Lt=L0;
+				else
+					Li= ( dipole_elements(j,m)*kron_delta(k,n)+dipole_elements(n,k)*kron_delta(j,m) )*E_field(w3i);
+					Lt=Li;
+				endif
+				Lt=-im_one/hbar*Lt*kron_delta(w_jk-w_iner,w_mn); % above if should be done only if kron_delta is not zero
+				Lt+=-im_one*w_jk*kron_delta(w_iner,w_jk);
+				L(p,s)=Lt;
+			endif
+		endfor
+	endfor
+endfor
+
+
+% now generally rhoL_dot=0=L*rhoL has infinite number of solutions
+% since we always can resclale rho vector with arbitrary constant
+% lets constrain our density matrix with some physical meaning
+% sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1
+% we will replace first row of Liouville operator with this condition
+% thus rhoLiouville_dot(1)=1
+for i=1:N
+	w2i=rhoLiouville_w(i);
+	m=rhoLiouville_r(i);
+	n=rhoLiouville_c(i);
+	w=modulation_freq(w2i);
+	if ((w==0) & (m==n))
+		L(1,i)=1;
+	endif
+endfor
+
+rhoLiouville_dot=rhoLiouville*0;
+% sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1
+% we will replace first row of Liouville operator with this condition
+% thus rhoLiouville_dot(1)=1
+rhoLiouville_dot(1)=1;
+
+
+%solving for density matrix vector
+rhoLiouville=L\rhoLiouville_dot;
+
+% this function create from Liouville density vector 
+% the density matrix with given modulation frequency
+function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels, Nfreq)
+	rho=zeros(Nlevels);
+	for r=1:Nlevels
+		for c=1:Nlevels
+			rho(r,c)=rhoLiouville((freqIndex-1)*Nlevels^2+(r-1)*Nlevels+c);
+		endfor
+	endfor
+endfunction
+
+
+rho_0=rhoOfFreq(rhoLiouville, 1, Nlevels, Nfreq)
+%rho_l=rhoOfFreq(rhoLiouville, Nfreq, Nlevels, Nfreq)
+
+%L*rhoLiouville
+