this code is for PSR calculation, added pseudo code to not recalculate atom properties matrices if they already saved

9 files changed, 967 insertions, 0 deletions

diff --git a/psr/basis_transformation.m b/psr/basis_transformation.m
new file mode 100644
index 0000000..82945fe
--- /dev/null
+++ b/psr/basis_transformation.m
@@ -0,0 +1,57 @@
+1;
+
+% matrix of circular to linear transformation
+% [x, y, z]' = lin2circ * [r, l, z]'
+function transformation_matrix = circ2lin()
+	transformation_matrix = ...
+		[ ...
+			[  1/sqrt(2),  1/sqrt(2), 0]; ...
+			[-1i/sqrt(2), 1i/sqrt(2), 0]; ...
+			[          0,          0, 1] ...
+		];
+endfunction
+
+
+% matrix of linear to circular transformation
+% [r, l, z]' = lin2circ * [x, y, z]'
+function transformation_matrix = lin2circ()
+	transformation_matrix = ...
+		[ ...
+			[  1/sqrt(2),  1i/sqrt(2), 0]; ...
+			[  1/sqrt(2), -1i/sqrt(2), 0]; ...
+			[          0,          0,  1] ...
+		];
+endfunction
+% linear basis rotation
+% x axis untouched
+% z and y rotated by angle theta around 'x' axis
+% [x_new, y_new, z_new]' = oldlin2newlin * [x_old, y_old, z_old]'
+function oldlin2newlin_m = oldlin2newlin(theta)
+	oldlin2newlin_m = [ ...
+			[ 1,          0,           0]; ...
+			[ 0, cos(theta), -sin(theta)]; ...
+			[ 0, sin(theta),  cos(theta)]...
+		];
+endfunction
+
+% rotate x polarized light by angle phi around
+% light propagation axis (Z)
+function [E_field_x, E_field_y] = rotXpolarization(phi, E_field_linear) 
+	% important negative frequency behave as they rotate in opposite direction
+	E_field_x=cos(phi)*E_field_linear;
+	E_field_y=sin(phi)*E_field_linear;
+endfunction					
+
+% transform  x,y,z linearly polarized light in the lab/light system coordinate
+% to left, right, linear along z atom system of coordinate
+% atom magnetic field is along new axis Z wich is at angle theta with respect to
+% light propagation direction
+function E_field_pos_freq=xyz_lin2atomic_axis_polarization(theta, E_field_lab_pos_freq)
+	coord_transf_m =  lin2circ() * oldlin2newlin(theta);
+	E_field_pos_freq.right  = coord_transf_m(1,1)*E_field_lab_pos_freq.x + coord_transf_m(1,2)*E_field_lab_pos_freq.y + coord_transf_m(1,3)*E_field_lab_pos_freq.z;
+	E_field_pos_freq.left   = coord_transf_m(2,1)*E_field_lab_pos_freq.x + coord_transf_m(2,2)*E_field_lab_pos_freq.y + coord_transf_m(2,3)*E_field_lab_pos_freq.z;
+	E_field_pos_freq.linear = coord_transf_m(3,1)*E_field_lab_pos_freq.x + coord_transf_m(3,2)*E_field_lab_pos_freq.y + coord_transf_m(3,3)*E_field_lab_pos_freq.z;
+endfunction
+	
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/psr/dipole_elementRb87D1line.m b/psr/dipole_elementRb87D1line.m
new file mode 100644
index 0000000..9235274
--- /dev/null
+++ b/psr/dipole_elementRb87D1line.m
@@ -0,0 +1,177 @@
+function d=dipole_elementRb87D1line(Fl,ml,Fu,mu)
+% Fl, ml  are F and m quantum numbers of lower state
+% Fu, mu  are F and m quantum numbers of upper state
+% F is total momentum and m is projection
+	d.left   = 0; %default return value
+	d.linear = 0; %default return value
+	d.right  = 0; %default return value
+	if ( mu==(ml+1) )
+		% sigma plus polarization
+		% ------ Fl=2 -> Fu=2 --------
+		if ( (ml==-2) & (Fl==2) & (Fu==2) )
+			d.right=sqrt(1/6);
+		endif
+		if ( (ml==-1) & (Fl==2) & (Fu==2) )
+			d.right=sqrt(1/4);
+		endif
+		if ( (ml== 0) & (Fl==2) & (Fu==2) )
+			d.right=sqrt(1/4);
+		endif
+		if ( (ml== 1) & (Fl==2) & (Fu==2) )
+			d.right=sqrt(1/6);
+		endif
+		if ( (ml== 2) & (Fl==2) & (Fu==2) )
+			d.right=0;
+		endif
+		% ------ Fl=2 -> Fu=1 --------
+		if ( (ml==-2) & (Fl==2) & (Fu==1) )
+			d.right=sqrt(1/2);
+		endif
+		if ( (ml==-1) & (Fl==2) & (Fu==1) )
+			d.right=sqrt(1/4);
+		endif
+		if ( (ml== 0) & (Fl==2) & (Fu==1) )
+			d.right=sqrt(1/12);
+		endif
+		if ( (ml== 1) & (Fl==2) & (Fu==1) )
+			d.right=0;
+		endif
+		if ( (ml== 2) & (Fl==2) & (Fu==1) )
+			d.right=0;
+		endif
+		% ------ Fl=1 -> Fu=2 --------
+		if ( (ml==-1) & (Fl==1) & (Fu==2) )
+			d.right = -sqrt(1/12);
+		endif
+		if ( (ml== 0) & (Fl==1) & (Fu==2) )
+			d.right = -sqrt(1/4);
+		endif
+		if ( (ml== 1) & (Fl==1) & (Fu==2) )
+			d.right = -sqrt(1/2);
+		endif
+		% ------ Fl=1 -> Fu=1 --------
+		if ( (ml==-1) & (Fl==1) & (Fu==1) )
+			d.right = -sqrt(1/12);
+		endif
+		if ( (ml== 0) & (Fl==1) & (Fu==1) )
+			d.right = -sqrt(1/12);
+		endif
+		if ( (ml== 1) & (Fl==1) & (Fu==1) )
+			d.right = 0;
+		endif
+	endif
+	if ( mu==(ml+0) )
+		% pi polarization
+		% ------ Fl=2 -> Fu=2 --------
+		if ( (ml==-2) & (Fl==2) & (Fu==2) )
+			d.linear=-sqrt(1/3);
+		endif
+		if ( (ml==-1) & (Fl==2) & (Fu==2) )
+			d.linear=-sqrt(1/12);
+		endif
+		if ( (ml== 0) & (Fl==2) & (Fu==2) )
+			d.linear=0;
+		endif
+		if ( (ml== 1) & (Fl==2) & (Fu==2) )
+			d.linear=sqrt(1/12);
+		endif
+		if ( (ml== 2) & (Fl==2) & (Fu==2) )
+			d.linear=sqrt(1/3);
+		endif
+		% ------ Fl=2 -> Fu=1 --------
+		if ( (ml==-2) & (Fl==2) & (Fu==1) )
+			d.linear = 0;
+		endif
+		if ( (ml==-1) & (Fl==2) & (Fu==1) )
+			d.linear=sqrt(1/4);
+		endif
+		if ( (ml== 0) & (Fl==2) & (Fu==1) )
+			d.linear=sqrt(1/3);
+		endif
+		if ( (ml== 1) & (Fl==2) & (Fu==1) )
+			d.linear=sqrt(1/4);
+		endif
+		if ( (ml== 2) & (Fl==2) & (Fu==1) )
+			d.linear = 0;
+		endif
+		% ------ Fl=1 -> Fu=2 --------
+		if ( (ml==-1) & (Fl==1) & (Fu==2) )
+			d.linear =  sqrt(1/4);
+		endif
+		if ( (ml== 0) & (Fl==1) & (Fu==2) )
+			d.linear =  sqrt(1/2);
+		endif
+		if ( (ml== 1) & (Fl==1) & (Fu==2) )
+			d.linear =  sqrt(1/4);
+		endif
+		% ------ Fl=1 -> Fu=1 --------
+		if ( (ml==-1) & (Fl==1) & (Fu==1) )
+			d.linear =  sqrt(1/12);
+		endif
+		if ( (ml== 0) & (Fl==1) & (Fu==1) )
+			d.linear = 0;
+		endif
+		if ( (ml== 1) & (Fl==1) & (Fu==1) )
+			d.linear = -sqrt(1/12);
+		endif
+	endif
+	if ( mu==(ml-1) )
+		% sigma minus polarization
+		% ------ Fl=2 -> Fu=2 --------
+		if ( (ml==-2) & (Fl==2) & (Fu==2) )
+			d.left = 0;
+		endif
+		if ( (ml==-1) & (Fl==2) & (Fu==2) )
+			d.left = -sqrt(1/6);
+		endif
+		if ( (ml== 0) & (Fl==2) & (Fu==2) )
+			d.left = -sqrt(1/4);
+		endif
+		if ( (ml== 1) & (Fl==2) & (Fu==2) )
+			d.left = -sqrt(1/4);
+		endif
+		if ( (ml== 2) & (Fl==2) & (Fu==2) )
+			d.left = -sqrt(1/6);
+		endif
+		% ------ Fl=2 -> Fu=1 --------
+		if ( (ml==-2) & (Fl==2) & (Fu==1) )
+			d.left = 0;
+		endif
+		if ( (ml==-1) & (Fl==2) & (Fu==1) )
+			d.left = 0;
+		endif
+		if ( (ml== 0) & (Fl==2) & (Fu==1) )
+			d.left =  sqrt(1/12);
+		endif
+		if ( (ml== 1) & (Fl==2) & (Fu==1) )
+			d.left =  sqrt(1/4);
+		endif
+		if ( (ml== 2) & (Fl==2) & (Fu==1) )
+			d.left =  sqrt(1/2);
+		endif
+		% ------ Fl=1 -> Fu=2 --------
+		if ( (ml==-1) & (Fl==1) & (Fu==2) )
+			d.left = -sqrt(1/2);
+		endif
+		if ( (ml== 0) & (Fl==1) & (Fu==2) )
+			d.left = -sqrt(1/4);
+		endif
+		if ( (ml== 1) & (Fl==1) & (Fu==2) )
+			d.left = -sqrt(1/12);
+		endif
+		% ------ Fl=1 -> Fu=1 --------
+		if ( (ml==-1) & (Fl==1) & (Fu==1) )
+			d.left = 0;
+		endif
+		if ( (ml== 0) & (Fl==1) & (Fu==1) )
+			d.left =  sqrt(1/12);
+		endif
+		if ( (ml== 1) & (Fl==1) & (Fu==1) )
+			d.left =  sqrt(1/12);
+		endif
+	endif
+endfunction
+
+
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/psr/field_description.m b/psr/field_description.m
new file mode 100644
index 0000000..738ebe3
--- /dev/null
+++ b/psr/field_description.m
@@ -0,0 +1,14 @@
+1;
+
+%EM field definition
+Ed=10.12; %drive
+Edc=conj(Ed);
+Ep=0.8*Ed; %probe
+Epc=conj(Ep);
+Em=-Ep; % opposite sideband (resulting from EOM modulation of drive)
+Emc=conj(Em);
+%wd=w13;
+%wp=w12;
+%wm=wd-(wp-wd);
+%light_positive_freq  = [wp, wd,  wp-wd];
+E_field               = [Ep, Ed,  0    ];
diff --git a/psr/output_results.m b/psr/output_results.m
new file mode 100644
index 0000000..eca6085
--- /dev/null
+++ b/psr/output_results.m
@@ -0,0 +1,31 @@
+1;
+
+
+load '/tmp/relative_transmission_vs_detuning.mat' ;
+
+figure(1);
+hold off;
+zoom_factor=1;
+plot(detuning_freq, zoom_factor*(relative_transmission_vs_detuning), '-');  
+title("relative transmission");
+xlabel("two photon detuning");
+
+%load 'xi_vs_detuning.mat' ;
+
+%figure(1); 
+	%hold off;
+	%plot(detuning_freq, imag(xi_linear), '-1;linear;');  
+	%hold on;
+	%plot(detuning_freq, imag(xi_left),   '-2;left;');  
+	%plot(detuning_freq, imag(xi_right),  '-3;right;');  
+	%title("probe absorption");
+	%hold off;
+%figure(2); 
+	%hold off;
+	%plot(detuning_freq, real(xi_linear), '-1;linear;');  
+	%hold on;
+	%plot(detuning_freq, real(xi_left),   '-2;left;');  
+	%plot(detuning_freq, real(xi_right),  '-3;right;');  
+	%title("probe dispersion");
+	%hold off;
+
diff --git a/psr/output_xi_results.m b/psr/output_xi_results.m
new file mode 100644
index 0000000..7dc3855
--- /dev/null
+++ b/psr/output_xi_results.m
@@ -0,0 +1,46 @@
+1;
+
+
+load '/tmp/xi_vs_detuning.mat' ;
+
+figure(1);
+hold off;
+plot(detuning_freq, real(xi_left-xi_right), '-');  
+title("differential real xi");
+xlabel("two photon detuning");
+
+figure(2);
+hold off;
+plot(detuning_freq, imag(xi_left-xi_right), '-');  
+title("differential imag xi");
+xlabel("two photon detuning");
+
+figure(3);
+hold off;
+plot(detuning_freq, imag(xi_left), '-', detuning_freq, imag(xi_right), '-');  
+title("imag xi");
+xlabel("two photon detuning");
+
+figure(4);
+hold off;
+plot(detuning_freq, real(xi_left), '-', detuning_freq, real(xi_right), '-');  
+title("real xi");
+xlabel("two photon detuning");
+
+%figure(1); 
+	%hold off;
+	%plot(detuning_freq, imag(xi_linear), '-1;linear;');  
+	%hold on;
+	%plot(detuning_freq, imag(xi_left),   '-2;left;');  
+	%plot(detuning_freq, imag(xi_right),  '-3;right;');  
+	%title("probe absorption");
+	%hold off;
+%figure(2); 
+	%hold off;
+	%plot(detuning_freq, real(xi_linear), '-1;linear;');  
+	%hold on;
+	%plot(detuning_freq, real(xi_left),   '-2;left;');  
+	%plot(detuning_freq, real(xi_right),  '-3;right;');  
+	%title("probe dispersion");
+	%hold off;
+
diff --git a/psr/psr.m b/psr/psr.m
new file mode 100644
index 0000000..7c8f604
--- /dev/null
+++ b/psr/psr.m
@@ -0,0 +1,149 @@
+1;
+clear all; 
+t0 = clock (); % we will use this latter to calculate elapsed time
+
+
+% load useful functions;
+useful_functions;
+
+% some physical constants
+useful_constants;
+
+basis_transformation; % load subroutines
+
+% load atom energy levels and decay description
+rb87_D1_line;
+%four_levels_with_polarization;
+%four_levels;
+%three_levels;
+%two_levels;
+
+% load EM field description
+field_description;
+
+%Nfreq=length(modulation_freq);
+
+
+
+%tune probe frequency
+detuning_p=0;
+N_detun_steps=100;
+%detuning_p_min=-B_field*gmg*4; % span  +/-4 Zeeman splitting
+detuning_p_min=-10040.0;
+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);
+detun_step=(detuning_p_max-detuning_p_min)/N_detun_steps;
+
+fprintf (stderr, "calculating atom properties\n");
+fflush (stderr);
+pfile='rb87_D1_line.m'; % the parent file from which L0_and_polarization_submatrices calculated
+cfile='L0m_and_polarizability_calculated.mat'; % the child  file to which calculated matrices writen 
+[s, err, msg] = stat (pfile); 
+if(err) 
+	%file does not exist
+	disp('Big troubles are coming, no file to define Hamiltonian)');
+	msg=cstrcat('File: ', pfile, ' is missing...exiting');
+	disp(msg);
+	return;
+else
+	pfile_mtime=s.mtime;
+endif
+[s, err, msg] = stat (cfile); 
+if(err) 
+	%file does not exist
+	cfile_mtime=0; 
+else 
+	cfile_mtime=s.mtime; 
+endif; 
+if ( cfile_mtime >= pfile_mtime)
+		% matrices already calculated and up to date, all we need to load them
+		load(cfile);
+	else
+		% calculate E_field independent properties of the 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 ...
+		);
+		save(cfile, 'L0m', 'polarizability_m');
+	endif
+elapsed_time = etime (clock (), t0);
+fprintf (stderr, "elapsed time so far is %.3f sec\n",elapsed_time);
+fflush (stderr);
+
+global atom_properties;
+atom_properties.L0m=L0m;
+atom_properties.polarizability_m=polarizability_m;
+atom_properties.dipole_elements=dipole_elements;
+
+%light_positive_freq  = [wp];
+E_field_drive         = [0    ];
+E_field_probe         = [Ep   ];
+E_field_zero          = [0    ];
+E_field_lab_pos_freq.linear = E_field_zero + (1.00000+0.00000i)*E_field_probe  + (1.00000+0.00000i)*E_field_drive; 
+%E_field_lab_pos_freq.right  = E_field_zero + (0.00000+0.00000i)*E_field_probe  + (0.00000+0.00000i)*E_field_drive;
+%E_field_lab_pos_freq.left   = E_field_zero + (0.00000+0.00000i)*E_field_probe  + (0.00000+0.00000i)*E_field_drive;
+
+% phi is angle between linear polarization and axis x
+phi=pi*2/8;
+% theta is angle between lab z axis (light propagation direction) and magnetic field axis (z')
+theta=0;
+% psi_el is the ellipticity parameter (phase difference between left and right polarization)
+psi_el=-3/180*pi;
+
+% we define light as linearly polarized 
+% where phi is angle between light polarization and axis x
+	% only sign of modulation frequency is important now
+	% we define actual frequency later on
+	[E_field_lab_pos_freq.x, E_field_lab_pos_freq.y] = rotXpolarization(phi, E_field_lab_pos_freq.linear);
+	E_field_lab_pos_freq.z=E_field_zero;
+
+	E_field_pos_freq=xyz_lin2atomic_axis_polarization(theta, E_field_lab_pos_freq);
+	E_field_pos_freq.left*=exp(I*psi_el);
+
+
+fprintf (stderr, "tuning laser in forloop to set conditions vs detuning\n");
+fflush (stderr);
+for detuning_p_cntr=1:N_detun_steps+1;
+	wp0=w_pf1-w_sf2;  %Fg=2 -> Fe=1
+	%wd=w_pf1-w_hpf_ground;
+	detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1);
+	wp=wp0+detuning_p;
+	light_positive_freq=[ wp];
+	% we calculate dc and negative frequiencies as well as amplitudes
+	[modulation_freq, E_field] = ...
+		light_positive_frequencies_and_amplitudes2full_set_of_modulation_frequencies_and_amlitudes(...
+			light_positive_freq, E_field_pos_freq);
+	freq_index=freq2index(wp,modulation_freq);
+
+	atom_field_problem.E_field          = E_field;
+	atom_field_problem.modulation_freq  = modulation_freq;
+	atom_field_problem.freq_index       = freq_index;
+
+	problems_cell_array{detuning_p_cntr}=atom_field_problem;
+
+
+	%kappa_p(detuning_p_cntr)=susceptibility_steady_state_at_freq( atom_field_problem);
+	detuning_freq(detuning_p_cntr)=detuning_p;
+endfor
+
+save '/tmp/problem_definition.mat' problems_cell_array atom_properties  detuning_freq ;
+fprintf (stderr, "now really hard calculations begin\n");
+fflush (stderr);
+% once we define all problems the main job is done here
+%kappa_p=cellfun( @susceptibility_steady_state_at_freq, problems_cell_array);
+%kappa_p=parcellfun(2, @susceptibility_steady_state_at_freq, problems_cell_array);
+[xi_linear, xi_left, xi_right]=parcellfun(2, @susceptibility_steady_state_at_freq, problems_cell_array);
+%relative_transmission_vs_detuning=parcellfun(2, @total_relative_transmission, problems_cell_array);
+%relative_transmission_vs_detuning=cellfun(@total_relative_transmission, problems_cell_array);
+
+save '/tmp/xi_vs_detuning.mat' detuning_freq xi_linear xi_left xi_right  ;
+%save '/tmp/relative_transmission_vs_detuning.mat' detuning_freq relative_transmission_vs_detuning;
+
+output_xi_results;
+
+elapsed_time = etime (clock (), t0)
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/psr/rb87_D1_line.m b/psr/rb87_D1_line.m
new file mode 100644
index 0000000..ecdd6a3
--- /dev/null
+++ b/psr/rb87_D1_line.m
@@ -0,0 +1,150 @@
+1;
+useful_constants;
+
+% note all frequency values are in MHz
+
+%  87Rb D1 line
+%
+%     m=-2   m=-1   m=0    m=1    m=2
+%     ----   ----   ----   ----   ----   |P,F=2>
+%
+%            ----   ----   ----          |P,F=1>
+%            m=-1   m=0    m=1        
+%
+%
+%
+%
+%
+%     m=-2   m=-1   m=0    m=1    m=2
+%     ----   ----   ----   ----   ----   |S,F=2>
+%
+%            ----   ----   ----          |S,F=1>
+%            m=-1   m=0    m=1        
+
+
+w_hpf_ground=6834;
+w_hpf_exited=817;
+w_sf2 = w_hpf_ground; %            Distance from |S,F=1> to |S,F=2>
+w_pf1       =1e9;  % something big Distance from |S,F=1> to |P,F=1>
+w_pf2       = w_pf1+w_hpf_exited; %Distance from |S,F=1> to |P,F=2> 
+gmg=.7; % gyro magnetic ration for ground level
+gme=.23; % gyro magnetic ration for exited level % CHECKME
+
+zeeman_splitting=0.0;
+B_field=zeeman_splitting/gmg;
+
+%bottom level   |F=1>
+levels( 1)=struct( "ang_momentum",  0, "total_momentum",   1, "m",  -1,  "energy",       0, "gm",  -gmg);
+levels( 2)=struct( "ang_momentum",  0, "total_momentum",   1, "m",   0,  "energy",       0, "gm",  -gmg);
+levels( 3)=struct( "ang_momentum",  0, "total_momentum",   1, "m",   1,  "energy",       0, "gm",  -gmg);
+
+%second bottom level |F=2>
+levels( 4)=struct( "ang_momentum",  0, "total_momentum",   2, "m",  -2,  "energy",   w_sf2, "gm",   gmg);
+levels( 5)=struct( "ang_momentum",  0, "total_momentum",   2, "m",  -1,  "energy",   w_sf2, "gm",   gmg);
+levels( 6)=struct( "ang_momentum",  0, "total_momentum",   2, "m",   0,  "energy",   w_sf2, "gm",   gmg);
+levels( 7)=struct( "ang_momentum",  0, "total_momentum",   2, "m",   1,  "energy",   w_sf2, "gm",   gmg);
+levels( 8)=struct( "ang_momentum",  0, "total_momentum",   2, "m",   2,  "energy",   w_sf2, "gm",   gmg);
+
+% first exited level |F=1>
+levels( 9)=struct( "ang_momentum",  1, "total_momentum",   1, "m",  -1,  "energy",   w_pf1, "gm",  -gme);
+levels(10)=struct( "ang_momentum",  1, "total_momentum",   1, "m",   0,  "energy",   w_pf1, "gm",  -gme);
+levels(11)=struct( "ang_momentum",  1, "total_momentum",   1, "m",   1,  "energy",   w_pf1, "gm",  -gme);
+
+% second exited level |F=2>
+levels(12)=struct( "ang_momentum",  1, "total_momentum",   2, "m",  -2,  "energy",   w_pf2, "gm",   gme);
+levels(13)=struct( "ang_momentum",  1, "total_momentum",   2, "m",  -1,  "energy",   w_pf2, "gm",   gme);
+levels(14)=struct( "ang_momentum",  1, "total_momentum",   2, "m",   0,  "energy",   w_pf2, "gm",   gme);
+levels(15)=struct( "ang_momentum",  1, "total_momentum",   2, "m",   1,  "energy",   w_pf2, "gm",   gme);
+levels(16)=struct( "ang_momentum",  1, "total_momentum",   2, "m",   2,  "energy",   w_pf2, "gm",   gme);
+
+Nlevels=size(levels)(2);
+
+
+H0=zeros(Nlevels);
+energy            = [1:Nlevels].*0;
+ang_momentum      = [1:Nlevels].*0;
+total_momentum    = [1:Nlevels].*0;
+m                 = [1:Nlevels].*0;
+gm                = [1:Nlevels].*0;
+
+for i=1:Nlevels
+	energy(i)         = levels(i).energy;
+	ang_momentum(i)   = levels(i).ang_momentum;
+	total_momentum(i) = levels(i).total_momentum;
+	m(i)              = levels(i).m;
+	gm(i)             = levels(i).gm;
+endfor
+H0=diag(energy)*hbar;
+
+
+dipole_elements.left   = zeros(Nlevels);
+dipole_elements.right  = zeros(Nlevels);
+dipole_elements.linear = zeros(Nlevels);
+% define dipole elements for transition from level j->k
+for j=1:Nlevels
+	for k=1:Nlevels
+		if ( abs(ang_momentum(j) - ang_momentum(k)) == 1)
+			%transition allowed for L =L' +/- 1
+			% but they correct within a factor, but not sign
+			if ( ( H0(j,j) < H0(k,k)) )
+				de=dipole_elementRb87D1line(...
+					total_momentum(j),m(j), ...
+					total_momentum(k),m(k)  ...
+				);
+				dipole_elements.left(j,k)   = de.left;
+				dipole_elements.right(j,k)  = de.right;
+				dipole_elements.linear(j,k) = de.linear;
+			endif	
+		endif
+	endfor
+endfor
+dipole_elements.linear+=dipole_elements.linear';
+dipole_elements.left+=dipole_elements.left';
+dipole_elements.right+=dipole_elements.right';
+
+maximum_dipole_elements=abs(dipole_elements.linear) + abs(dipole_elements.left) + abs(dipole_elements.right);
+
+%defasing matrix
+g_deph=0;
+g_dephasing=zeros(Nlevels);
+% dephasing only for non zero dipole elements (am I right?)
+g_dephasing=g_deph*(abs(maximum_dipole_elements) != 0);
+
+% decay matrix g(i,j) correspnds to decay from i-->j
+gamma=6;
+g_decay=zeros(Nlevels);
+for i=1:Nlevels
+	for j=1:Nlevels
+		if (  H0(i,i) > H0(j,j) )
+			% only upper levels decaying and decay is always positive
+			g_decay(i,j)=gamma * abs( maximum_dipole_elements(i,j) );
+		endif
+	endfor
+endfor
+
+% ground level mixing need to be artificial
+gamma_hpf=.0001;
+for i=1:Nlevels
+	for j=1:Nlevels
+		% it would be better to introduce a level corresponding to the bath
+		% but this slows calculations
+		% So 
+		% ground hyperfine are equally mixed together
+		if (  (abs( H0(i,i) - H0(j,j)) == w_hpf_ground*hbar ) || ...
+			    (abs( H0(i,i) - H0(j,j)) == 0 ) )
+			% i and j correspond to 2 ground levels
+			% we cannot decay to itself
+			if ( i != j )	
+				% total eight ground levels F=2 and F=1 so we decay in 7 channels
+				g_decay(i,j)+=gamma_hpf/7; 
+			endif
+		endif
+	endfor
+endfor
+
+% apply B field Zeeman splitting
+energy+=B_field * ( gm .* m);
+% convert frequency to energy units
+H0=diag(energy)*hbar;
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/psr/useful_constants.m b/psr/useful_constants.m
new file mode 100644
index 0000000..82d1d1a
--- /dev/null
+++ b/psr/useful_constants.m
@@ -0,0 +1,6 @@
+1;
+
+% some physical constants
+hbar=1;
+im_one=0+1i;
+el_charge=1;
diff --git a/psr/useful_functions.m b/psr/useful_functions.m
new file mode 100644
index 0000000..eb3e880
--- /dev/null
+++ b/psr/useful_functions.m
@@ -0,0 +1,337 @@
+1;
+
+function ret=decay_total(g_decay,i)
+% calculate total decay for particular level taking in account all branches
+	ret=sum(g_decay(i,:));
+endfunction
+
+function ret=kron_delta(i,j)
+% kroneker delta symbol
+	if ((i==j))
+		ret=1;
+	else
+		ret=0;
+	endif
+endfunction
+
+function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels)
+% this function create from Liouville density vector 
+% the density matrix with given modulation frequency
+	rho=zeros(Nlevels);
+	rho(:)=rhoLiouville((freqIndex-1)*Nlevels^2+1:(freqIndex)*Nlevels^2);
+	rho=rho.';
+endfunction
+
+function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq)
+% 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)
+	N  = Nfreq*Nlevels*Nlevels;
+	rho_size = Nlevels*Nlevels;
+	rhoLiouville_w=zeros(N,1);
+	rhoLiouville_r=zeros(N,1);
+	rhoLiouville_c=zeros(N,1);
+	
+	w=1:Nfreq;
+	w_tmplate=(repmat(w,rho_size,1))(:);
+	rhoLiouville_w=w_tmplate;
+	r=1:Nlevels;
+	r_tmplate=(repmat(r,Nlevels,1))(:);
+	rhoLiouville_r=(repmat(r_tmplate,Nfreq,1))(:)';
+	c=(1:Nlevels)';% hold column value of rho_rc
+	rhoLiouville_c=repmat(c,Nfreq*Nlevels,1);
+endfunction	
+
+
+function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
+		Nlevels, ...
+		H0, g_decay, g_dephasing, dipole_elements ...
+		)
+% create (Nlevels*Nlevels)x*(Nlevels*Nlevels)
+% sub matrices of Liouville operator 
+% which repeat themselves for each modulation frequency
+% based on recipe from Eugeniy Mikhailov thesis
+	%-------------------------
+	useful_constants;
+	rho_size=Nlevels*Nlevels;
+
+	% now we create Liouville indexes list
+	[Ndummy, rhoLiouville_w_notused, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,1);
+
+	kron_delta_m=eye(Nlevels);
+	% note that L0 and decay parts depend only on combination of indexes
+	% jk,mn but repeats itself for every frequency
+	L0m=zeros(rho_size); % (Nlevels^2)x(Nlevels^2) matrix
+	decay_part_m=zeros(rho_size); % (NxN)x(NxN) matrix
+	% polarization matrix will be multiplied by field amplitude letter
+	% polarization is part of perturbation part of Hamiltonian
+	polarizability_m.linear = zeros(rho_size); % (NxN)x(NxN) matrix
+	polarizability_m.left   = zeros(rho_size); % (NxN)x(NxN) matrix
+	polarizability_m.right  = zeros(rho_size); % (NxN)x(NxN) matrix
+	for p=1:rho_size
+		% p= j*Nlevels+k 
+		% this might speed up stuff since less matrix passed back and force
+		j=rhoLiouville_r(p);
+		k=rhoLiouville_c(p);
+		for s=1:rho_size
+			% s= m*Nlevels+n 
+			m=rhoLiouville_r(s);
+			n=rhoLiouville_c(s);
+
+			% calculate unperturbed part (Hamiltonian without EM field)
+			L0m(p,s)=H0(j,m)*kron_delta_m(k,n)-H0(n,k)*kron_delta_m(j,m);
+			decay_part_m(p,s)= ...
+				( ...
+					  decay_total(g_decay,k)/2 ...
+					+ decay_total(g_decay,j)/2 ...
+					+ g_dephasing(j,k) ...
+				)* kron_delta_m(j,m)*kron_delta_m(k,n) ...
+				- kron_delta_m(m,n)*kron_delta_m(j,k)*g_decay(m,j) ;
+			polarizability_m.linear(p,s)= ( dipole_elements.linear(j,m)*kron_delta_m(k,n)-dipole_elements.linear(n,k)*kron_delta_m(j,m) );
+			polarizability_m.left(p,s)= ( dipole_elements.left(j,m)*kron_delta_m(k,n)-dipole_elements.left(n,k)*kron_delta_m(j,m) );
+			polarizability_m.right(p,s)= ( dipole_elements.right(j,m)*kron_delta_m(k,n)-dipole_elements.right(n,k)*kron_delta_m(j,m) );
+		endfor
+	endfor
+	L0m=-im_one/hbar*L0m - decay_part_m; 
+endfunction
+
+function L=Liouville_operator_matrix( ...
+		N, ... 
+		L0m, polarizability_m, ...
+		E_field, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
+		)
+% Liouville operator matrix construction
+% based on recipe from Eugeniy Mikhailov thesis
+	%-------------------------
+	useful_constants;
+	L=zeros(N); % NxN matrix
+	Nfreq=length(modulation_freq);
+
+	% Lets be supper smart and speed up L matrix construction
+	% since 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;
+	
+	% creating building blocks of L by rho_size * rho_size
+	for w3i=1:Nfreq
+		w_iner=modulation_freq(w3i);
+		if ((w_iner == 0))
+			% calculate unperturbed part (Hamiltonian without EM field)
+			L_sub{w3i}=L0m;
+		else
+			% calculate perturbed part (Hamiltonian with EM field)
+			% in other word interactive part of Hamiltonian 
+			L_sub{w3i} = ...
+				-im_one/hbar*polarizability_m.linear * E_field.linear(w3i) ...
+				-im_one/hbar*polarizability_m.left   * E_field.left(w3i)   ...
+				-im_one/hbar*polarizability_m.right  * E_field.right(w3i)  ...
+				;
+		endif
+	endfor
+
+	% Liouville matrix operator has Nlevels*Nlevels blocks 
+	% which governed by the same modulation frequency
+	for p_freq_cntr=1:Nfreq
+		p0=1+(p_freq_cntr-1)*rho_size;
+		% we guaranteed to know frequency of final and initial rhoLiouville
+		w1i=rhoLiouville_w(p0); % final
+		w_jk=modulation_freq(w1i);
+		for s_freq_cntr=1:Nfreq
+			s0=1+(s_freq_cntr-1)*rho_size;
+			w2i=rhoLiouville_w(s0); % initial
+			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
+			w3i=(w_l == modulation_freq);
+			if (any(w3i))	
+				% yey, requested modulation frequency exist
+				% lets do  L sub matrix filling
+				% at most we should have only one matching frequency
+				w_iner=modulation_freq(w3i);
+				L(p0:p0+rho_size-1,s0:s0+rho_size-1) = L_sub{w3i};
+
+			endif
+		endfor
+		% diagonal elements are self modulated
+		% due to rotating wave approximation
+		L(p0:p0+rho_size-1,p0:p0+rho_size-1)+= -im_one*w_jk*eye(rho_size);
+	endfor
+endfunction
+
+
+function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ...
+		N, L, ...
+		modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
+% 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= 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
+
+function kappa=susceptibility(wi, rhoLiouville, dipole_elements, E_field)
+% calculate susceptibility for the field at given frequency index
+	Nlevels=( size(dipole_elements.linear)(1) );
+	rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
+	kappa.linear=0;
+	kappa.left=0;
+	kappa.right=0;
+
+	kappa.linear = sum( sum( transpose(dipole_elements.linear) .* rho ) );
+	kappa.left   = sum( sum( transpose(dipole_elements.left)   .* rho ) );
+	kappa.right  = sum( sum( transpose(dipole_elements.right)  .* rho ) );
+
+	kappa.linear /= E_field.linear( wi ); 
+	kappa.left   /= E_field.left( wi );
+	kappa.right  /= E_field.right( wi ); 
+endfunction
+
+function index=freq2index(freq, modulation_freq) 
+% convert modulation freq to its index in the modulation_freq vector
+	index=[1:length(modulation_freq)](modulation_freq==freq);
+endfunction
+
+function rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq)
+% calculates rhoLiouville vector assuming steady state situation and normalization of rho_ii to 1
+	Nlevels=sqrt( size(L0m)(1) );
+	Nfreq=length(modulation_freq);
+
+	% now we create Liouville indexes list
+	[N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
+
+	% 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;
+endfunction
+
+function [xi_linear, xi_left, xi_right]=susceptibility_steady_state_at_freq( atom_field_problem)
+	% find steady state susceptibility at particular modulation frequency element
+	% at given E_field 
+	global atom_properties;
+	L0m                = atom_properties.L0m              ; 
+	polarizability_m   = atom_properties.polarizability_m ; 
+	dipole_elements    = atom_properties.dipole_elements  ; 
+
+	E_field            = atom_field_problem.E_field          ; 
+	modulation_freq    = atom_field_problem.modulation_freq  ; 
+	freq_index         = atom_field_problem.freq_index       ; 
+
+	rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
+	xi=susceptibility(freq_index, rhoLiouville, dipole_elements, E_field);
+	xi_linear = xi.linear;
+	xi_right  = xi.right;
+	xi_left   = xi.left;
+endfunction		
+
+function [dEdz_coef_linear, dEdz_coef_left, dEdz_coef_right] =dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field)
+% complex absorption coefficient 
+% at given E_field frequency component
+	Nlevels=( size(dipole_elements.linear)(1) );
+	rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
+
+	dEdz_coef_linear = sum( sum( transpose(dipole_elements.linear) .* rho ) );
+	dEdz_coef_left   = sum( sum( transpose(dipole_elements.left)   .* rho ) );
+	dEdz_coef_right  = sum( sum( transpose(dipole_elements.right)  .* rho ) );
+
+endfunction 
+
+function [dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem)
+% complex absorption coefficient for each field
+	global atom_properties;
+	L0m                = atom_properties.L0m              ; 
+	polarizability_m   = atom_properties.polarizability_m ; 
+	dipole_elements    = atom_properties.dipole_elements  ; 
+
+	E_field            = atom_field_problem.E_field          ; 
+	modulation_freq    = atom_field_problem.modulation_freq  ; 
+	Nfreq=length(modulation_freq);
+
+	rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
+	for wi=1:Nfreq
+	 [dEdz_linear(wi), dEdz_left(wi), dEdz_right(wi)] = dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field);
+	endfor
+endfunction		
+
+function relative_transmission=total_relative_transmission(atom_field_problem)
+% summed across all frequencies field absorption coefficient
+	global atom_properties;
+	[dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem);
+	%dEdz_total= dEdz_linear .* conj(dEdz_linear) ...
+									 %+dEdz_left   .* conj(dEdz_left) ...
+									 %+dEdz_right  .* conj(dEdz_right);
+	%total_absorption = sum(dEdz_total);
+
+	%total_absorption = |E|^2 - | E - dEdz | ^2		
+	E_field.linear = atom_field_problem.E_field.linear;
+	E_field.right  = atom_field_problem.E_field.right;
+	E_field.left   = atom_field_problem.E_field.left;
+	
+	modulation_freq    = atom_field_problem.modulation_freq  ; 
+	pos_freq_indexes=(modulation_freq>0);
+	% WARNING INTRODUCED HERE BUT MUST BE DEFINE IN ATOM PROPERTIES FILE
+	%n_atoms - proportional to a number of interacting atoms
+	n_atoms=500;
+	transmited_intensities_vector = ...	
+		 abs(E_field.linear + n_atoms*(1i)*dEdz_linear).^2 ...
+		+abs(E_field.right  + n_atoms*(1i)*dEdz_right).^2 ...
+		+abs(E_field.left   + n_atoms*(1i)*dEdz_left).^2;
+	transmited_intensities_vector = transmited_intensities_vector(pos_freq_indexes);
+	input_intensities_vector =  ...	
+		 abs(E_field.linear).^2   ...
+		+abs(E_field.right).^2    ...
+		+abs(E_field.left).^2;
+	input_intensities_vector = input_intensities_vector(pos_freq_indexes);
+		
+	relative_transmission = sum(transmited_intensities_vector) / sum(input_intensities_vector);
+endfunction
+
+% create full list of atom modulation frequencies from positive light frequency amplitudes
+function [modulation_freq, E_field_amplitudes] = ...
+light_positive_frequencies_and_amplitudes2full_set_of_modulation_frequencies_and_amlitudes(...
+	positive_light_frequencies, positive_light_field_amplitudes )
+	% we should add 0 frequency as first element of our frequencies list
+	modulation_freq    = cat(2, 0, positive_light_frequencies, -positive_light_frequencies);
+	% negative frequencies have complex conjugated light fields amplitudes
+	E_field_amplitudes.left = cat(2, 0, positive_light_field_amplitudes.left, conj(positive_light_field_amplitudes.left) );
+	E_field_amplitudes.right = cat(2, 0, positive_light_field_amplitudes.right, conj(positive_light_field_amplitudes.right) );
+	E_field_amplitudes.linear = cat(2, 0, positive_light_field_amplitudes.linear, conj(positive_light_field_amplitudes.linear) );
+endfunction	
+
+
+
+% vim: ts=2:sw=2:fdm=indent