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