1;
clear all;

t0 = clock (); % we will use this latter to calculate elapsed time


% load useful functions;
useful_functions;

% some physical constants
useful_constants;

% load atom energy levels and decay description
%four_levels;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nlevels=4;
w1=1e9;
w2=0;
w_hpf=6800;
w3=w_hpf;
w4=w_hpf+.1; % separation of levels |3> and |4> somewhat like Zeeman splitting
w12=w1-w2;
w13=w1-w3;

%
%      -----------  |1>
%         /     \
%    E_d /       \ 
%       /         \ E_p
%      /           \
%    -------- |3>   \
%    --------   |4>  \
%                     \
%               ___________ |2>


% unperturbed Hamiltonian energy levels
levels_energy=[ w1, 0,  w3, w4];
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_23=.001;
g_decay=zeros(Nlevels);
g_decay(1,2)=gamma; %upper level decay
g_decay(1,3)=gamma; %upper level decay
g_decay(1,4)=gamma; %upper level decay
g_decay(3,2)=gamma_23; % lower levels mixing
g_decay(2,3)=gamma_23; % lower levels mixing
g_decay(4,2)=gamma_23; % lower levels mixing
g_decay(2,4)=gamma_23; % 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);
g_dephasing(1,4)=g_deph;
g_dephasing(4,1)=g_dephasing(4,1);



% dipole matrix
dipole_elements=zeros(Nlevels);
dipole_elements(1,2)=1;
dipole_elements(2,1)=dipole_elements(1,2);
dipole_elements(1,3)=1;
dipole_elements(3,1)=dipole_elements(1,3);
dipole_elements(1,4)=1;
dipole_elements(4,1)=dipole_elements(1,4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%EM field definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ep=0.01; %probe
Epc=conj(Ep);
Ed=.1; %drive
Edc=conj(Ed);
Em=-Ep; % opposite sideband (resulting from EOM modulation of drive)
Emc=conj(Em);
%wd=w13;
%wp=w12;
%wm=wd-(wp-wd);
%modulation_freq=[0,  wp, wd, wm,   -wp, -wd, -wm,  wp-wd, wd-wp];
E_field        =[0,  Ep, Ed, Em,   Epc, Edc, Emc,  0,     0    ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Nfreq=length(modulation_freq);



%tune probe frequency
detuning_p=0;
N_detun_steps=100;
detuning_p_min=-1;
detuning_p_max=-detuning_p_min;
detuning_freq=zeros(1,N_detun_steps+1);
kappa_p      =zeros(1,N_detun_steps+1);
kappa_m      =zeros(1,N_detun_steps+1);

w_pf1=1e9;
w_hpf_ground=6800;
wd=w_pf1-w_hpf_ground;
wp0=w_pf1;
wp=wp0+detuning_p_min;
wm=wd-(wp-wd);
%modulation_freq=[0,  wp, wd, wm,   -wp, -wd, -wm,  wp-wd, wd-wp];
%E_field        =[0,  Ep, Ed, Em,   Epc, Edc, Emc,  0,     0    ];
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
[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
rhoLiouville=zeros(N,1);

dipoles=dipole_elements;
Efld=E_field;
clear dipole_elements;
dipole_elements.right=dipoles;
dipole_elements.left=dipoles;
dipole_elements.linear=dipoles;
clear E_field;
E_field.right=Efld;
E_field.left=Efld;
E_field.linear=Efld;
% calculate E_field independent properties of athe atom
% to be used as sub matrix templates for Liouville operator matrix
[L0m, polarizability_m]=L0_and_polarization_submatrices( ...
		Nlevels, ...
		H0, g_decay, g_dephasing, dipole_elements ...
		);
% Liouville operator matrix construction
L=Liouville_operator_matrix( 
		N, 
		L0m, polarizability_m,
		E_field, 
		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
		);


%use the fact that sum(rho_ii)=1 to constrain solution
[rhoLiouville_dot, L]=constrain_rho_and_match_L(
		N, L,
		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c);


%solving for density matrix vector
rhoLiouville=L\rhoLiouville_dot;


L_new=L;
rhoLiouville_new=rhoLiouville;

% uncomment to update reference file
%save 'L_and_rhoL_referenced.mat' L rhoLiouville

clear L rhoLiouville;

load 'L_and_rhoL_referenced.mat' 
diff_with_reference_L=sum(sum(abs(L_new-L)))

diff_with_reference_rhoL=(sum(abs(rhoLiouville_new -rhoLiouville)))

elapsed_time = etime (clock (), t0)