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);

	% Lets be supper smart and speed up L matrix construction
	% it has a lot of voids.
	% By creation of rhoLiouville we know that there are
	% consequent chunks of rho_ij modulated with same frequency
	% this means that rhoLiouville is split in to Nfreq chunks 
	% with length Nlevels*Nlevels=N/Nfreq
	rho_size=N/Nfreq;
	for p_freq_cntr=1:Nfreq
		for s_freq_cntr=1:Nfreq
			p0=1+(p_freq_cntr-1)*rho_size;
			s0=1+(s_freq_cntr-1)*rho_size;
			% we guaranteed to know frequency of final and initial rhoLiouville
			w1i=rhoLiouville_w(p0);
			w2i=rhoLiouville_w(s0);
			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 frequency in the list of available frequencies
			% but since we not guaranteed to find it lets assign temporary 0 to Liouville matrix element
			decay_part=0;
			Lt=0;
			for w3i=1:Nfreq
				% at most we should have only one matching frequency
				w_iner=modulation_freq(w3i);
				if ((w_iner == w_l))
					% yey, requested modulation frequency exist
					% lets do  L submatrix filling
					for p_cntr=0:(rho_size-1)
						p=p0+p_cntr;
						for s_cntr=0:(rho_size-1)
							s=s0+s_cntr;

							decay_part=0;
							Lt=0;
							
							j=rhoLiouville_r(p);
							k=rhoLiouville_c(p);
							m=rhoLiouville_r(s);
							n=rhoLiouville_c(s);

							%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; 
							% diagonal elements are self modulated
							if ((p == s))
								Lt+=-im_one*w_jk;
							endif
							Lps=Lt;
							L(p,s)=Lps;
						endfor
					endfor
				endif
			endfor
		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

% calculate sucseptibility for the field at given frequency index
function kappa=sucseptibility(wi, rhoLiouville, dipole_elements, Nlevels, Nfreq)

	rho=rhoOfFreq(rhoLiouville, wi, Nlevels, Nfreq);
	kappa=0;
	for i=1:Nlevels
		for j=1:Nlevels
			kappa+=dipole_elements(j,i)*rho(i,j);
		endfor
	endfor
endfunction