% some physical constants
hbar=1;
im_one=0+1i;

%three_levels;
two_levels;

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)
N=length(modulation_freq)*Nlevels*Nlevels;
rhoLiouville=zeros(N,1);
rhoLiouville_w=rhoLiouville;
rhoLiouville_r=rhoLiouville;
rhoLiouville_c=rhoLiouville;
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


% 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)
	if ((i==j))
		ret=1;
	else
		ret=0;
	endif
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;
		decay_part=0;
		Lt=0;
		for w3i=1:Nfreq
			w_iner=modulation_freq(w3i);
			decay_part=0;
			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) ;
					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
				% no need for above kron_delta since the same conditon checked in the outer if statement
				Lt=-im_one/hbar*Lt - decay_part; 
			endif
		endfor
		if ((p == s))
			Lt+=-im_one*w_jk;
		endif
		L(p,s)=Lt;
	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;
	else
		L(1,i)=0;
	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_1=rhoOfFreq(rhoLiouville, 2, Nlevels, Nfreq)
rho_2=rhoOfFreq(rhoLiouville, 3, Nlevels, Nfreq)
%rho_l=rhoOfFreq(rhoLiouville, Nfreq, Nlevels, Nfreq)

%L*rhoLiouville