1;

% calculate total decay for particular level taking in account all branches
function ret=decay_total(g_decay,i)
	ret=0;
	for k=1:size(g_decay)(1)
		ret=ret+g_decay(i,k);
	endfor
endfunction

% kroneker delta symbol
function ret=kron_delta(i,j)
	if ((i==j))
		ret=1;
	else
		ret=0;
	endif
endfunction

% 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

% we unwrap density matrix to Liouville density vector and assign all possible
% modulation frequencies as well
% resulting vector should be Nlevels x Nlevels x length(modulation_freq)
function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq)
	N=Nfreq*Nlevels*Nlevels;
	rhoLiouville_w=zeros(N,1);
	rhoLiouville_r=zeros(N,1);
	rhoLiouville_c=zeros(N,1);
	i=0;
	for w=1:Nfreq;
		for r=1:Nlevels
			for c=1:Nlevels
				i+=1;
				rhoLiouville(i)=0;
				rhoLiouville_w(i)=w; % hold frequency modulation index
				rhoLiouville_r(i)=r; % hold row value of rho_rc
				rhoLiouville_c(i)=c; % hold column value of rho_rc
			endfor
		endfor
	endfor
endfunction	


% Liouville operator matrix construction
function L=Liouville_operator_matrix( 
		N, 
		H0, g_decay, g_dephasing, dipole_elements,
		E_field, 
		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
		)
	%-------------------------
	useful_constants;
	L=zeros(N); % NxN matrix

	Nfreq=length(modulation_freq);

	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))
						% calculate unperturbed part (Hamiltonian without EM field)
						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
						% calculate perturbed part (Hamiltonian with EM field)
						% in othe word interactive part of Hamiltonian 
						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
endfunction


% 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
function [rhoLiouville_dot, L]=constran_rho_and_match_L(
		N, L,
		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
	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= zeros(N,1);
	% 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;
endfunction