Merge branch 'master' of qo.physics.wm.edu:Nresonances

32 files changed, 161 insertions, 1014 deletions

diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm/Makefile
index 578f2b0..f8b3559 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm/Makefile
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/Makefile
@@ -15,7 +15,7 @@ M_FILES = $(patsubst %.xmds,%.m,$(XMDS_FILES))
 XMDS = xmds2
 XSIL2GRAPHICS = xsil2graphics
 
-all: $(M_FILES) 
+all: $(RUN_FILES)
 
 %.run: %.xmds
 	$(XMDS) $<
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/Nlevels_no_dopler_with_z_4wm.xmds b/xmds2/Nlevels_no_dopler_with_z_4wm/Nlevels_no_dopler_with_z_4wm.xmds
index 44b13f7..bb0feb2 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm/Nlevels_no_dopler_with_z_4wm.xmds
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/Nlevels_no_dopler_with_z_4wm.xmds
@@ -36,7 +36,7 @@
 
 		VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to
 		normalize it to something else look drho/dt equation.
-		No need to renormalizes eta as long as its express through i
+		No need to renormalizes eta as long as its express through
 		the upper level decay rate in the same units as Rabi frequency.
 	</description>
 
@@ -49,7 +49,8 @@
 				const double lambda=794.7e-9; //wavelength in m
 				const double Kvec = 2*M_PI/lambda;  // k-vector
 				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
-				const double eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				// eta will be calculated in the <arguments> section
+				double eta = 0;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
 				
 				//  ---------  Atom and cell properties -------------------------
 				// range of Maxwell distribution atomic velocities
@@ -57,9 +58,10 @@
 				// above mass expression is written as (expression is isotopic_mass * atomic_mass_unit)
 
 				// Average sqrt(v^2) in Maxwell distribution for one dimension
-				const double v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell related parameters will be calculated in <arguments> section
+				double v_thermal_averaged=0;
 				// Maxwell distribution velocities range to take in account in [m/s]
-				const double V_maxwell_min = -4*v_thermal_averaged, V_maxwell_max = -V_maxwell_min; // there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
+				double V_maxwell_min = 0, V_maxwell_max = 0;
 
 				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
 				const double gt=0.01 *(2*M_PI*1e6);
@@ -101,6 +103,18 @@
 					 I am guessing detunings are too large and thus it became a stiff equation-->
 			<!--! make sure it is not equal to zero!-->
 			<argument name="Temperature"  type="real" default_value="5" />
+			<!-- This will be executed after arguments/parameters are parsed -->
+			<!-- Read the code Luke: took me a while of reading the xmds2 sources to find it -->
+			<![CDATA[
+				// Average sqrt(v^2) in Maxwell distribution for one dimension
+				v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell distribution velocities range to take in account in [m/s]
+				// there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
+				V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; 
+
+				// eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;
+			]]>
     </arguments>
 		<bing />
 		<fftw plan="patient" />
@@ -112,7 +126,7 @@
 	<geometry>
 		<propagation_dimension> z </propagation_dimension>
 		<transverse_dimensions>
-			<dimension name="t"   lattice="1000"   domain="(-14.0e-7, 4.0e-7)" />
+			<dimension name="t"   lattice="1000"   domain="(-0.2e-6, 0.4e-6)" />
 		</transverse_dimensions>
 	</geometry>
 
@@ -310,7 +324,7 @@
 		</group>
 
 		<group>
-      <sampling basis="t(100)" initial_sample="yes">
+      <sampling basis="t(1000)" initial_sample="yes">
 				<dependencies>density_matrix</dependencies>
 				<moments>
 					r11_out r22_out r33_out r44_out 
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/fast_light/plot_fields_propagation_I4.gp b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_light/plot_fields_propagation_I4.gp
new file mode 120000
index 0000000..134a380
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_light/plot_fields_propagation_I4.gp
@@ -0,0 +1 @@
+../plot_fields_propagation_I4.gp
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/Makefile
index 697fd27..fcc817c 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/Makefile
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/Makefile
@@ -9,11 +9,17 @@ GNUPLOT_FILES = $(wildcard *.gp)
 XSIL2GRAPHICS = xsil2graphics
 
 # fast light
-# PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=3.8e7 --E4o=6.3e4
+PARAMS =  \
+	  --Ndens=1e15 \
+	  --Lcell=10.0e-2 \
+	  --Temperature=1e-9 \
+	  --Pwidth=0.1e-6 \
+	  --delta1=0 --delta2=0 --delta3=0 \
+	  --E1o=1.9e7 --E2o=1e2 --E3o=3.8e7 --E4o=0
+
+
 # slow light EIT
 #PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=0 --E4o=0
-#Fast light to Slow light switch
-PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=2e7 --E2o=3e3 --E3o=6e6 --E4o=3e3
 
 all: $(XSIL_FILES) Nlevels_no_dopler_with_z_4wm.xsil $(M_FILES) plot png
 
@@ -23,9 +29,12 @@ Nlevels_no_dopler_with_z_4wm.xsil: ../Nlevels_no_dopler_with_z_4wm.run
 %.m: %.xsil
 	$(XSIL2GRAPHICS) $<
 
-plot: $(M_FILES) $(GNUPLOT_FILES) 
-	octave pp_I2.m
+pretty_plots: plot $(GNUPLOT_FILES)
 	gnuplot plot_fields_propagation_I2.gp
+	gnuplot plot_fields_propagation_I4.gp
+
+plot: $(M_FILES)
+	octave pp_I2.m
 
 clean:
 	rm -f $(CC_FILES) $(RUN_FILES) $(M_FILES) $(XSIL_FILES) *.wisdom.fftw3 *.dat octave-core *.wisdom *.pdf
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/README b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/README
index 70f7262..ec7d200 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/README
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/README
@@ -1,5 +1,10 @@
 fast to slow transition
-PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=3.8e7 --E4o=6.3e4
+PARAMS =  \
+	  --Ndens=1e15 \
+	  --Lcell=10.0e-2 \
+	  --Temperature=1e-9 \
+	  --delta1=0 --delta2=0 --delta3=0 \
+	  --E1o=1.9e7 --E2o=3.1e5 --E3o=3.8e7 --E4o=0 
 
 1.5 cm seems to give largest advance
 2.5 cm gives zero delay
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/plot_fields_propagation_I4.gp b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/plot_fields_propagation_I4.gp
new file mode 120000
index 0000000..134a380
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/plot_fields_propagation_I4.gp
@@ -0,0 +1 @@
+../plot_fields_propagation_I4.gp
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/pp_I2.m
index 73487e0..4a1ed9b 100644..120000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/pp_I2.m
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/fast_to_slow_switch/pp_I2.m
@@ -1,67 +1 @@
-Nlevels_no_dopler_with_z_4wm
-
-%% field propagation
-z_1=z_1*100; % z in cm
-t_1=t_1*1e6; % time now measured in uS
-figure(1)
-set(gca,'fontsize',20); 
-imagesc(z_1, t_1, I2_out_1); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_2')
-title('I_2')
-
-xskip=1;
-yskip=10;
-map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
-
-
-
-print('-color','fields_propagation_I2.eps')
-
-
-
-%% fields before and after the cell
-figure(2)
-%set(gca,'fontsize',30); 
-plot( ...
-	t_1,I2_out_1(:,1),'.-;before;', "linewidth", 4,  ...
-	t_1,I2_out_1(:,end), '-;after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I_2 (1/s)^2')
-title('I_2 before and after cell')
-legend('location', 'southwest');
-
-[b, a]=butter(3, 0.05);
-I2_out_after=I2_out_1(:,end);
-I2_out_after_filtered=filtfilt(b,a,I2_out_after);
-settling_time=0.8; %uS
-t_good_indx=t_1> min(t_1 + settling_time);
-[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
-delay_time=t_1(max_pos_after)-t_1(max_pos_before);
-printf('Second field delay time  = %f uS\n',delay_time);
-
-%set(gca,'fontsize',40); 
-set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','fields_before_after_cell_I2.eps')
-
-figure(4)
-I2_max_in=max(I2_out_1(t_good_indx,1));
-I2_max_out=max(I2_out_1(t_good_indx,end));
-I2_in_norm=(I2_out_1(:,1))/I2_max_in;
-I2_out_norm=(I2_out_1(:,end))/I2_max_out;
-tmin=-0.05;
-tmax=0.05;
-indx=(t_1>=tmin & t_1<=tmax);  % soom in in time to this region
-plot( ...
-	t_1(indx),I2_in_norm(indx),'.-;before;', "linewidth", 4,  ...
-	t_1(indx),I2_out_norm(indx), '-;after;', "linewidth", 4 ...
-	)
-xlim([tmin,tmax],'manual');	
-xlabel('t (uS)')
-ylabel('I_2')
-title('I_2 before and after cell normalized')
-set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','probe_before_after_cell_I2_normalized.eps')
-
+../pp_I2.m
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/plot_fields_propagation_I2.gp b/xmds2/Nlevels_no_dopler_with_z_4wm/plot_fields_propagation_I2.gp
index b6721d4..77b4870 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm/plot_fields_propagation_I2.gp
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/plot_fields_propagation_I2.gp
@@ -12,4 +12,4 @@ set ylabel "t ({/Symbol m}S)"
 set zlabel "I_2 (1/S)"
 set nokey
 #set view map
-splot [0:1.5][-0.4:0.4] 'I2.dat' 
+splot [0:][-0.2:0.2] 'I2.dat' 
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/plot_fields_propagation_I2.gp b/xmds2/Nlevels_no_dopler_with_z_4wm/plot_fields_propagation_I4.gp
index 77b4870..3849e2f 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/plot_fields_propagation_I2.gp
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/plot_fields_propagation_I4.gp
@@ -1,5 +1,5 @@
 set terminal postscript portrait enhanced color solid size 5,3.5
-set output 'fields_propagation_I2.eps'
+set output 'fields_propagation_I4.eps'
 set dgrid3d 100,100 qnorm 4
 set pm3d map
 #set contour
@@ -12,4 +12,4 @@ set ylabel "t ({/Symbol m}S)"
 set zlabel "I_2 (1/S)"
 set nokey
 #set view map
-splot [0:][-0.2:0.2] 'I2.dat' 
+splot [0:][-0.2:0.2] 'I4.dat' 
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm/pp_I2.m
index 45d4913..3e5c174 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm/pp_I2.m
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/pp_I2.m
@@ -1,6 +1,6 @@
-Nlevels_no_dopler_with_z_4wm
+Nlevels_no_dopler_with_z_4wm;
 
-%% field propagation
+%% field I2  propagation
 z_1=z_1*100; % z in cm
 t_1=t_1*1e6; % time now measured in uS
 figure(1)
@@ -13,33 +13,50 @@ xlabel('z (cm)')
 ylabel('t (uS)')
 zlabel('I_2')
 title('I_2')
+print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I2.eps')
 
-xskip=1;
-yskip=10;
-%map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
+desired_x_size=200;
+desired_y_size=200;
+xskip=ceil(length(z_1)/desired_x_size);
+yskip=ceil(length(t_1)/desired_y_size);
+map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
 
 
+%% field I4  propagation
+figure(2)
+
+imagesc(z_1, t_1, I4_out_1); colorbar
+tmin=-0.4;
+tmax= 0.4;
+ylim([tmin,tmax],'manual');	
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_4')
+title('I_4')
+print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I4.eps')
+map2dat('I4.dat',z_1,t_1, I4_out_1, xskip, yskip);
+
 
-print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I2.eps')
 
 
 
 %% fields before and after the cell
-figure(2)
+figure(3)
 %set(gca,'fontsize',30); 
 plot( ...
-	t_1,I2_out_1(:,1),'.-;before;', "linewidth", 4,  ...
-	t_1,I2_out_1(:,end), '-;after;', "linewidth", 4 ...
+	t_1,I2_out_1(:,1),'.-;I_2 before;', "linewidth", 4  ...
+	,t_1,I2_out_1(:,end), '-;I_2 after;', "linewidth", 4 ...
+	,t_1,I4_out_1(:,end), '-;I_4 after;', "linewidth", 4 ...
 	)
 xlabel('t (uS)')
-ylabel('I_2 (1/s)^2')
-title('I_2 before and after cell')
-legend('location', 'southwest');
+ylabel('I (1/s)^2')
+title('Fields before and after cell')
+legend('location', 'northeast');
 
 [b, a]=butter(3, 0.05);
 I2_out_after=I2_out_1(:,end);
 I2_out_after_filtered=filtfilt(b,a,I2_out_after);
-settling_time=0.8; %uS
+settling_time=0.01; %uS
 t_good_indx=t_1> min(t_1 + settling_time);
 [m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
 delay_time=t_1(max_pos_after)-t_1(max_pos_before);
@@ -47,25 +64,29 @@ printf('Second field delay time  = %f uS\n',delay_time);
 
 %set(gca,'fontsize',40); 
 %set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell_I2.eps')
+print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell.eps')
 
 figure(4)
 I2_max_in=max(I2_out_1(t_good_indx,1));
 I2_max_out=max(I2_out_1(t_good_indx,end));
+I4_max_out=max(I4_out_1(t_good_indx,end));
 I2_in_norm=(I2_out_1(:,1))/I2_max_in;
 I2_out_norm=(I2_out_1(:,end))/I2_max_out;
-tmin=-0.05;
-tmax=0.05;
+I4_out_norm=(I4_out_1(:,end))/I4_max_out;
+tmin=-.05;
+tmax=.05;
 indx=(t_1>=tmin & t_1<=tmax);  % soom in in time to this region
 plot( ...
-	t_1(indx),I2_in_norm(indx),'.-;before;', "linewidth", 4,  ...
-	t_1(indx),I2_out_norm(indx), '-;after;', "linewidth", 4 ...
+	t_1(indx),I2_in_norm(indx),'.-;I_2 before;', "linewidth", 4,  ...
+	t_1(indx),I2_out_norm(indx), '-;I_2 after;', "linewidth", 4 ...
+	,t_1(indx),I4_out_norm(indx), '-;I_4 after;', "linewidth", 4 ...
 	)
 legend('location', 'southeast');
 xlim([tmin,tmax],'manual');	
+ylim([0.5,1],'manual');	
 xlabel('t (uS)')
-ylabel('I_2')
-title('I_2 before and after cell normalized')
+ylabel('I normalized')
+title('Fields before and after cell normalized')
 %set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','-depsc2', '-tight','-S200,120', 'probe_before_after_cell_I2_normalized.eps')
+print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell_normalized.eps')
 
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/Makefile
index 95ebb0a..95ebb0a 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/Makefile
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/Makefile
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/map2dat.m b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/map2dat.m
index 14fae30..14fae30 120000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/map2dat.m
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/map2dat.m
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/plot_fields_propagation_I2.gp b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/plot_fields_propagation_I2.gp
new file mode 120000
index 0000000..e11c569
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/plot_fields_propagation_I2.gp
@@ -0,0 +1 @@
+../plot_fields_propagation_I2.gp
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/plot_fields_propagation_I4.gp b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/plot_fields_propagation_I4.gp
new file mode 120000
index 0000000..134a380
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/plot_fields_propagation_I4.gp
@@ -0,0 +1 @@
+../plot_fields_propagation_I4.gp
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/pp_I2.m
new file mode 120000
index 0000000..4a1ed9b
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/pulse_split/pp_I2.m
@@ -0,0 +1 @@
+../pp_I2.m
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/slow_light/plot_fields_propagation_I4.gp b/xmds2/Nlevels_no_dopler_with_z_4wm/slow_light/plot_fields_propagation_I4.gp
new file mode 120000
index 0000000..134a380
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/slow_light/plot_fields_propagation_I4.gp
@@ -0,0 +1 @@
+../plot_fields_propagation_I4.gp
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Makefile
deleted file mode 100644
index 578f2b0..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Makefile
+++ /dev/null
@@ -1,54 +0,0 @@
-###  -*- make -*-
-### This file is part of the Debian xmds package
-### Copyright (C) 2006 Rafael Laboissiere
-### This file is relased under the GNU General Public License
-### NO WARRANTIES!
-
-### This makefile can be used to build and run the XMDS examples
-
-XMDS_FILES = $(shell ls *.xmds)
-RUN_FILES = $(patsubst %.xmds,%.run,$(XMDS_FILES))
-CC_FILES = $(patsubst %.xmds,%.cc,$(XMDS_FILES))
-XSIL_FILES = $(patsubst %.xmds,%.xsil,$(XMDS_FILES))
-M_FILES = $(patsubst %.xmds,%.m,$(XMDS_FILES))
-
-XMDS = xmds2
-XSIL2GRAPHICS = xsil2graphics
-
-all: $(M_FILES) 
-
-%.run: %.xmds
-	$(XMDS) $<
-	mv $(patsubst %.xmds,%,$<) $@
-
-%.xsil: %.run
-	./$<
-
-%.m: %.xsil
-	$(XSIL2GRAPHICS) $<
-
-plot: $(M_FILES)
-	octave pp.m
-
-clean:
-	rm -f $(CC_FILES) $(RUN_FILES) $(M_FILES) $(XSIL_FILES) *.wisdom.fftw3 *.dat octave-core *.wisdom *.pdf 
-	rm -f $(png_targets)
-
-eps_targets  = $(wildcard *.eps)
-pdf_targets  = $(eps_targets:%.eps=%.pdf)	
-png_targets  = $(pdf_targets:%.pdf=%.png)
-
-png: pdf $(png_targets) 
-
-$(png_targets): %.png : %.pdf
-	convert -density 300 $< $@
-
-pdf: $(pdf_targets) 
-
-$(pdf_targets): %.pdf : %.eps
-	cat $< | ps2eps -B > __tt.eps
-	epspdf __tt.eps $@
-	rm -f __tt.eps
-	#ps2eps -B $< | epspdf $< $@
-.PRECIOUS: %.run %.xsil %.m
-.PHONY: all clean
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Nlevels_no_dopler_with_z_4wm.xmds b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Nlevels_no_dopler_with_z_4wm.xmds
deleted file mode 100644
index 5102c46..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Nlevels_no_dopler_with_z_4wm.xmds
+++ /dev/null
@@ -1,349 +0,0 @@
-<?xml version="1.0"?>
-<simulation xmds-version="2">
-
-	<name>Nlevels_no_dopler_with_z_4wm</name>
-
-	<author>Eugeniy Mikhailov</author>
-	<description>
-		License GPL.
-
-		Solving 4 level atom in N-field configuration,
-		with field propagation along spatial axis Z
-		no Doppler broadening
-
-		We assume four-wave mixing condition when w3-w4=w2-w1 i.e. fields E3 and E4 drive the same
-		resonance as fields E2 and E1.
-
-
-		*  --------------- |4>
-		*    \        \    
-		*     \ E3     \    -------- |3>  
-		*      \     E4 \   /    \
-		*       \        \ / E2   \ 
-		*        \        /        \ E1
-		*   |2> --------------      \
-		*                   \        \
-		*                    \        \
-		*                   ------------- |1>
-		*
-
-
-		We are solving 
-			dE/dz+(1/c)*dE/dt=i*eta*rho_ij,  where j level is higher then i.
-			Note that E is actually a Rabi frequency of electromagnetic field not the EM field
-		in xmds terms it looks like
-			dE_dz = i*eta*rhoij - 1/c*L[E], here we moved t dependence to Fourier space
-
-		VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to
-		normalize it to something else look drho/dt equation.
-		No need to renormalizes eta as long as its express through i
-		the upper level decay rate in the same units as Rabi frequency.
-	</description>
-
-	<features>
-		<globals>
-			<![CDATA[
-				const double pi = M_PI; 
-				const double c=3.e8;
-				const double k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K]
-				const double lambda=794.7e-9; //wavelength in m
-				const double Kvec = 2*M_PI/lambda;  // k-vector
-				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
-				const double eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
-				
-				//  ---------  Atom and cell properties -------------------------
-				// range of Maxwell distribution atomic velocities
-				const double mass = (86.909180527 * 1.660538921e-27); // atom mass in [kg] 
-				// above mass expression is written as (expression is isotopic_mass * atomic_mass_unit)
-
-				// Average sqrt(v^2) in Maxwell distribution for one dimension
-				const double v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
-				// Maxwell distribution velocities range to take in account in [m/s]
-				const double V_maxwell_min = -4*v_thermal_averaged, V_maxwell_max = -V_maxwell_min; // there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
-
-				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
-				const double gt=0.01 *(2*M_PI*1e6);
-				// Natural linewidth  of j's level in [1/s]
-				const double G3=3.0 *(2*M_PI*1e6);
-				const double G4=3.0 *(2*M_PI*1e6);
-
-				// total decay of i-th level branching ratios. Rij branching of i-th level to j-th
-				const double R41=0.5, R42=0.5;
-				const double R31=0.5, R32=0.5;
-
-
-				complex  E1c, E2c, E3c, E4c; // Complex conjugated Rabi frequencies
-
-				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
-			]]>
-		</globals>
-		<validation kind="run-time"/> <!--allows to put ranges as variables-->
-		<benchmark />
-    <arguments>
-			<!-- Rabi frequency divided by 2 in [1/s] -->
-      <argument name="E1o" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
-      <argument name="E2o" type="real" default_value="0.05*(2*M_PI*1e6)" />
-      <argument name="E3o" type="real" default_value="2*3.0*(2*M_PI*1e6)" />
-      <argument name="E4o" type="real" default_value=".01*(2*M_PI*1e6)" />
-			<!-- Fields detuning in [1/s] -->
-      <argument name="delta1"  type="real" default_value="0.0" />
-      <argument name="delta2"  type="real" default_value="0.0" />
-      <argument name="delta3"  type="real" default_value="0.0" />
-			<!--Pulse duration/width [s] -->
-			<argument name="Pwidth"  type="real" default_value="0.1e-6" />
-			<!--  Atom and cell properties -->
-			<!--Cell length [m] -->
-			<argument name="Lcell"  type="real" default_value="1.5e-2" />
-			<!--Density of atoms [1/m^3] -->
-			<argument name="Ndens"  type="real" default_value="1e15" />
-			<!--Atoms temperature [K] -->
-			<!--TODO: looks like Temperature > 10 K knocks solver, 
-					 I am guessing detunings are too large and thus it became a stiff equation-->
-			<!--! make sure it is not equal to zero!-->
-			<argument name="Temperature"  type="real" default_value="5" />
-    </arguments>
-		<bing />
-		<fftw plan="patient" />
-		<openmp />
-		<auto_vectorise />
-	</features>
-
-	<!-- 'z' and 't' to have dimensions [m] and [s]   -->
-	<geometry>
-		<propagation_dimension> z </propagation_dimension>
-		<transverse_dimensions>
-			<dimension name="t"   lattice="1000"   domain="(-0.2e-6, 0.4e-6)" />
-		</transverse_dimensions>
-	</geometry>
-
-	<!-- Rabi frequency --> 
-	<vector name="E_field" type="complex" initial_space="t">
-		<components>E1 E2 E3 E4</components>
-		<initialisation>
-			<![CDATA[
-			// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
-			// as well as time
-			E1=E1o;
-			E2=E2o*exp(-pow( ((t-0.0)/Pwidth),2) );
-			E3=E3o;
-			E4=E4o;
-			]]>
-		</initialisation>
-	</vector>
-
-	<vector name="density_matrix" type="complex" initial_space="t">
-		<!--<components>r11 r22 r33 r44 r12 r13 r14 r23 r24 r34   r21 r31 r41 r32 r42 r43</components>-->
-		<!--<components>r11 r22 r33 r44 r12 r13 r14 r23 r24 r34</components>-->
-		<components>r11 r22 r33 r12 r13 r14 r23 r24 r34 r44</components>
-		<!--
-				 note one of the level population is redundant since
-				 r11+r22+r33+r44=1
-				 so r11 is missing
-		-->
-		<initialisation>
-			<!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
-			<!-- Comment out no light field initial conditions 
-			<![CDATA[
-			// Note: 
-			// convergence is really slow if all populations concentrated at the bottom level |1>
-			// this is because if r11=1, everything else is 0 and then every small increment 
-			// seems to be huge and adaptive solver makes smaller and smaller steps.
-			// As quick and dirty fix I reshuffle initial population  
-			// so some of the population sits at the  second ground level |2>
-			// TODO: Fix above. Make the equation of motion for r11 
-			//       and express other level, let's say r44
-			//       through population normalization
-			r11 = 1; 
-			          r22 = 0; r33 = 0; r44 = 0;
-			 r12 = 0; r13 = 0; r14 = 0;
-			 r23 = 0; r24 = 0;
-			 r34 = 0;
-			]]>
-			-->
-			<!-- Below initialization assumes strong E1 and E3 which were shining for long time before
-					 we even start to look at the problem -->
-			<!-- Precalculated by Nate via Mathematica steady state solver -->
-				<dependencies>E_field</dependencies>
-				<![CDATA[
-						E1c = conj(E1);
-						E2c = conj(E2);
-						E3c = conj(E3);
-						E4c = conj(E4);
-
-				// IMPORTANT: assumes no detunings
-				r11 = (gt*(4*mod2(E1) + (G3 + gt)*(G3 + 2*gt))*((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt)) + 4*mod2(E3)*G4*(R41 -
-				R42)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r13 = -((E1c*gt*(G3 + gt)*i*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt))
-				+ 4*mod2(E3)*G4*(R41 -
-				R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r22 = (gt*(4*mod2(E3) + (G4 + gt)*(G4 + 2*gt))*((G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt)) + 4*mod2(E1)*G3*(-R31 +
-				R32)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r24 = -((E3c*gt*(G4 + gt)*i*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt))
-				+ 4*mod2(E1)*G3*(-R31 +
-				R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r33 = (2*mod2(E1)*gt*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt)) +
-				4*mod2(E3)*G4*(R41 -
-				R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42));
-
-				r44 = (2*mod2(E3)*gt*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt)) +
-				4*mod2(E1)*G3*(-R31 +
-				R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42));
-
-				]]>
-		</initialisation>
-	</vector>
-
-	<sequence>
-		<!--For this set of conditions ARK45 is faster than ARK89-->
-		<!--ARK45 is good for small detuning when all frequency like term are close to zero-->
-		<integrate algorithm="ARK45" tolerance="1e-5" interval="Lcell"> 
-		<!--RK4 is good for large detunings when frequency like term are big, it does not try to be too smart about adaptive step which ARK seems to make too small-->
-		<!--When ARK45 works it about 3 times faster then RK4 with 1000 steps-->
-		<!--<integrate algorithm="RK4" steps="100" interval="1.5e-2">-->
-		<!--SIC algorithm seems to be much slower and needs fine 'z'  step tuning and much finer time grid-->
-		<!--For example I had to quadruple the time grid from 1000 to 4000 when increased z distance from 0.02 to 0.04-->
-
-		<!--<integrate algorithm="SIC" interval="4e-2" steps="200">-->
-			<samples>200 200</samples>
-			<operators>
-        <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
-					<integration_vectors>density_matrix</integration_vectors>
-          <dependencies>E_field</dependencies>
-          <boundary_condition kind="left">
-						<!--This set boundary condition at all 'z'  and left border of 't' (i.e. min(t))-->
-						<!--
-            <![CDATA[
-							r11 = 0; r22 = 1; r33 = 0; r44 = 0;
-							r12 = 0; r13 = 0; r14 = 0;
-							r23 = 0; r24 = 0;
-							r34 = 0;
-							printf("z= %g, t= %g\n", z, t);
-            ]]>
-						-->
-          </boundary_condition>
-					<![CDATA[
-						E1c = conj(E1);
-						E2c = conj(E2);
-						E3c = conj(E3);
-						E4c = conj(E4);
-
-						r21=conj(r12);
-						r31=conj(r13);
-						r41=conj(r14);
-						r32=conj(r23);
-						r42=conj(r24);
-						r43=conj(r34);
-						//r44=1- r11 - r22 - r33;
-
-						// Equations of motions according to Nate's mathematica code
-						dr11_dt = gt/2 - gt*r11 + i*(-(E1*r13) - E4*r14 + E1c*r31 + E4c*r41) + G3*r33*R31 + G4*r44*R41;
-						dr12_dt = -(gt*r12) + i*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1c*r32 + E4c*r42);
-						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1c*r11) - E2c*r12 - delta1*r13 + E1c*r33 + E4c*r43);
-						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4c*r11) - E3c*r12 - (delta1 - delta2 + delta3)*r14 + E1c*r34 + E4c*r44);
-						dr22_dt = gt/2 - gt*r22 + i*(-(E2*r23) - E3*r24 + E2c*r32 + E3c*r42) + G3*r33*R32 + G4*r44*R42;
-						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1c*r21) - E2c*r22 - delta2*r23 + E2c*r33 + E3c*r43);
-						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4c*r21) - E3c*r22 - delta3*r24 + E2c*r34 + E3c*r44);
-						dr33_dt = i*(E1*r13 + E2*r23 - E1c*r31 - E2c*r32) - (G3 + gt)*r33;
-
-						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1*r14 + E2*r24 - E4c*r31 - E3c*r32 + (delta2 - delta3)*r34);
-						dr44_dt = i*(E4*r14 + E3*r24 - E4c*r41 - E3c*r42) - (G4 + gt)*r44;
-					]]>
-        </operator>
-				<operator kind="ex" constant="yes">
-					<operator_names>Lt</operator_names>
-					<![CDATA[
-					Lt = i*1./c*kt;
-					]]>
-				</operator>
-        <integration_vectors>E_field</integration_vectors>
-				<dependencies>density_matrix</dependencies>
-          <![CDATA[
-					dE1_dz = i*eta*conj(r13) -Lt[E1] ;
-					dE2_dz = i*eta*conj(r23) -Lt[E2] ;
-					dE3_dz = i*eta*conj(r24) -Lt[E3] ;
-					dE4_dz = i*eta*conj(r14) -Lt[E4] ;
-          ]]>
-			</operators>
-		</integrate>
-	</sequence>
-
-
-
-
-	<!-- The output to generate -->
-	<output format="binary" filename="Nlevels_no_dopler_with_z_4wm.xsil">
-		<group>
-      <sampling basis="t(1000)" initial_sample="yes">
-				<dependencies>E_field</dependencies>
-				<moments>I1_out I2_out I3_out I4_out</moments>
-				<![CDATA[
-				I1_out = mod2(E1);
-				I2_out = mod2(E2);
-				I3_out = mod2(E3);
-				I4_out = mod2(E4);
-				]]>
-			</sampling>
-		</group>
-
-		<group>
-      <sampling basis="t(1000)" initial_sample="yes">
-				<dependencies>density_matrix</dependencies>
-				<moments>
-					r11_out r22_out r33_out r44_out 
-					r12_re_out r12_im_out r13_re_out r13_im_out r14_re_out r14_im_out
-					                      r23_re_out r23_im_out r24_re_out r24_im_out 
-					                                            r34_re_out r34_im_out
-				</moments>
-				<![CDATA[
-				// populations output 
-				r11_out = r11.Re();
-				r22_out = r22.Re();
-				r33_out = r33.Re();
-				r44_out = r44.Re();
-				// coherences output 
-				r12_re_out = r12.Re();
-				r12_im_out = r12.Im();
-				r13_re_out = r13.Re();
-				r13_im_out = r13.Im();
-				r14_re_out = r14.Re();
-				r14_im_out = r14.Im();
-				r23_re_out = r23.Re();
-				r23_im_out = r23.Im();
-				r24_re_out = r24.Re();
-				r24_im_out = r24.Im();
-				r34_re_out = r34.Re();
-				r34_im_out = r34.Im();
-				]]>
-			</sampling>
-		</group>
-	</output>
-
-</simulation>
-	
-<!-- 
-vim: ts=2 sw=2 foldmethod=indent: 
--->
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/Makefile
deleted file mode 100644
index 7581c04..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/Makefile
+++ /dev/null
@@ -1,54 +0,0 @@
-###  -*- make -*-
-### This makefile can be used to build and run the XMDS examples
-
-
-XSIL_FILES = Nlevels_no_dopler_with_z_4wm.xsil
-M_FILES = $(patsubst %.xsil,%.m,$(XSIL_FILES))
-GNUPLOT_FILES = $(wildcard *.gp)
-
-XSIL2GRAPHICS = xsil2graphics
-
-# fast light
-PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=3.8e7 --E4o=0
-# slow light EIT
-#PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=0 --E4o=0
-
-all: $(XSIL_FILES) Nlevels_no_dopler_with_z_4wm.xsil $(M_FILES) plot png
-
-Nlevels_no_dopler_with_z_4wm.xsil: ../Nlevels_no_dopler_with_z_4wm.run
-	$< $(PARAMS) | grep "Time elapsed for simulation is:" > exact_analysis_execution_time.txt
-
-%.m: %.xsil
-	$(XSIL2GRAPHICS) $<
-
-plot: $(M_FILES) $(GNUPLOT_FILES) 
-	octave pp_I2.m
-	gnuplot plot_fields_propagation_I2.gp
-
-clean:
-	rm -f $(CC_FILES) $(RUN_FILES) $(M_FILES) $(XSIL_FILES) *.wisdom.fftw3 *.dat octave-core *.wisdom *.pdf
-	rm -f $(pdf_targets)
-	rm -f $(eps_targets)
-
-real_clean: clean
-	rm -f $(png_targets)
-
-eps_targets  = $(wildcard *.eps)
-pdf_targets  = $(eps_targets:%.eps=%.pdf)	
-png_targets  = $(pdf_targets:%.pdf=%.png)
-
-pdf: $(pdf_targets) 
-
-$(pdf_targets): %.pdf : %.eps
-	cat $< | ps2eps -B > __tt.eps
-	epspdf __tt.eps $@
-	rm -f __tt.eps
-	#ps2eps -B $< | epspdf $< $@
-
-png: pdf $(png_targets) 
-
-$(png_targets): %.png : %.pdf
-	convert -density 300 $< $@
-
-.PRECIOUS: %.run %.xsil %.m
-.PHONY: all clean
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/pp_I2.m
deleted file mode 100644
index c4b6333..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_to_slow_switch/pp_I2.m
+++ /dev/null
@@ -1,91 +0,0 @@
-Nlevels_no_dopler_with_z_4wm;
-
-%% field I2  propagation
-z_1=z_1*100; % z in cm
-t_1=t_1*1e6; % time now measured in uS
-figure(1)
-%set(gca,'fontsize',20); 
-imagesc(z_1, t_1, I2_out_1); colorbar
-tmin=-0.4;
-tmax= 0.4;
-ylim([tmin,tmax],'manual');	
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_2')
-title('I_2')
-print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I2.eps')
-
-desired_x_size=200;
-desired_y_size=200;
-xskip=ceil(length(z_1)/desired_x_size);
-yskip=ceil(length(t_1)/desired_y_size);
-map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
-
-
-%% field I4  propagation
-figure(2)
-
-imagesc(z_1, t_1, I4_out_1); colorbar
-tmin=-0.4;
-tmax= 0.4;
-ylim([tmin,tmax],'manual');	
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_2')
-title('I_2')
-print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I4.eps')
-
-
-
-
-
-%% fields before and after the cell
-figure(3)
-%set(gca,'fontsize',30); 
-plot( ...
-	t_1,I2_out_1(:,1),'.-;I_2 before;', "linewidth", 4  ...
-	,t_1,I2_out_1(:,end), '-;I_2 after;', "linewidth", 4 ...
-	,t_1,I4_out_1(:,end), '-;I_4 after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I (1/s)^2')
-title('Fields before and after cell')
-legend('location', 'northeast');
-
-[b, a]=butter(3, 0.05);
-I2_out_after=I2_out_1(:,end);
-I2_out_after_filtered=filtfilt(b,a,I2_out_after);
-settling_time=0.01; %uS
-t_good_indx=t_1> min(t_1 + settling_time);
-[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
-delay_time=t_1(max_pos_after)-t_1(max_pos_before);
-printf('Second field delay time  = %f uS\n',delay_time);
-
-%set(gca,'fontsize',40); 
-%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell.eps')
-
-figure(4)
-I2_max_in=max(I2_out_1(t_good_indx,1));
-I2_max_out=max(I2_out_1(t_good_indx,end));
-I4_max_out=max(I4_out_1(t_good_indx,end));
-I2_in_norm=(I2_out_1(:,1))/I2_max_in;
-I2_out_norm=(I2_out_1(:,end))/I2_max_out;
-I4_out_norm=(I4_out_1(:,end))/I4_max_out;
-tmin=-.05;
-tmax=.05;
-indx=(t_1>=tmin & t_1<=tmax);  % soom in in time to this region
-plot( ...
-	t_1(indx),I2_in_norm(indx),'.-;I_2 before;', "linewidth", 4,  ...
-	t_1(indx),I2_out_norm(indx), '-;I_2 after;', "linewidth", 4 ...
-	,t_1(indx),I4_out_norm(indx), '-;I_4 after;', "linewidth", 4 ...
-	)
-legend('location', 'southeast');
-xlim([tmin,tmax],'manual');	
-ylim([0.5,1],'manual');	
-xlabel('t (uS)')
-ylabel('I normalized')
-title('Fields before and after cell normalized')
-%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell_normalized.eps')
-
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/map2dat.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/map2dat.m
deleted file mode 120000
index 90fc33d..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/map2dat.m
+++ /dev/null
@@ -1 +0,0 @@
-../../matlab_helper_files/map2dat.m
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp.m
deleted file mode 100644
index b88fa6e..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp.m
+++ /dev/null
@@ -1,190 +0,0 @@
-Nlevels_no_dopler_with_z_4wm
-
-%% field propagation
-z_1=z_1*100; % z in cm
-t_1=t_1*1e6; % time now measured in uS
-figure(1)
-subplot(2,2,1); imagesc(z_1, t_1, I1_out_1); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_1')
-title('I_1')
-subplot(2,2,2); imagesc(z_1, t_1, I2_out_1); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_2')
-title('I_2')
-subplot(2,2,3); imagesc(z_1, t_1, I3_out_1); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_3')
-title('I_3')
-subplot(2,2,4); imagesc(z_1, t_1, I4_out_1); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_4')
-title('I_4')
-
-
-print('-color','fields_propagation.eps')
-
-
-
-%% fields before and after the cell
-figure(2)
-subplot(2,2,1);
-plot( ...
-	t_1,I1_out_1(:,1),'-;before;', "linewidth", 4,  ...
-	t_1,I1_out_1(:,end), '-;after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I_1 (1/s)^2')
-title('I_1 before and after cell')
-subplot(2,2,2);
-plot( ...
-	t_1,I2_out_1(:,1),'-;before;', "linewidth", 4,  ...
-	t_1,I2_out_1(:,end), '-;after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I_2 (1/s)^2')
-title('I_2 before and after cell')
-subplot(2,2,3);
-plot( ...
-	t_1,I3_out_1(:,1),'-;before;', "linewidth", 4,  ...
-	t_1,I3_out_1(:,end), '-;after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I_3 (1/s)^2')
-title('I_3 before and after cell')
-
-[b, a]=butter(3, 0.05);
-I2_out_after=I2_out_1(:,end);
-I2_out_after_filtered=filtfilt(b,a,I2_out_after);
-settling_time=0.8; %uS
-t_good_indx=t_1> min(t_1 + settling_time);
-[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
-delay_time=t_1(max_pos_after)-t_1(max_pos_before);
-printf('Second field delay time  = %f uS\n',delay_time);
-
-print('-color','fields_before_after_cell.eps')
-
-subplot(2,2,4);
-plot( ...
-	t_1,I4_out_1(:,1),'-;before;', "linewidth", 4,  ...
-	t_1,I4_out_1(:,end), '-;after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I_3 (1/s)^2')
-title('I_3 before and after cell')
-
-figure(4)
-I2_max_in=max(I2_out_1(t_good_indx,1));
-I2_max_out=max(I2_out_1(t_good_indx,end));
-I2_in_norm=(I2_out_1(:,1))/I2_max_in;
-I2_out_norm=(I2_out_1(:,end))/I2_max_out;
-tmin=-0.05;
-tmax=0.05;
-indx=(t_1>=tmin & t_1<=tmax);  % soom in in time to this region
-plot( ...
-	t_1(indx),I2_in_norm(indx),'-;before;', "linewidth", 4,  ...
-	t_1(indx),I2_out_norm(indx), '-;after;', "linewidth", 4 ...
-	)
-xlim([tmin,tmax],'manual');	
-xlabel('t (uS)')
-ylabel('I_2')
-title('I_2 before and after cell normalized')
-print('-color','probe_before_after_cell.eps')
-
-return;
-
-%% all density matrix elements in one plot
-% diagonal populations, 
-% upper triangle real part of coherences, 
-% lower diagonal imaginary part of coherences
-z_2=z_2*100; % z in cm
-t_2=t_2*1e6; % time now measured in uS
-figure(3)
-subplot(4,4,1); imagesc (z_2, t_2, r11_out_2);  caxis([0,1]); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('rho_{11}')
-title('rho_{11}')
-subplot(4,4,6); imagesc (z_2, t_2, r22_out_2);  caxis([0,1]); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('rho_{22}')
-title('rho_{22}')
-subplot(4,4,11); imagesc (z_2, t_2, r33_out_2);  caxis([0,1]); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('rho_{33}')
-title('rho_{33}')
-subplot(4,4,16); imagesc (z_2, t_2, r44_out_2); caxis([0,1]);  colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('rho_{44}')
-title('rho_{44}')
-% real parts of coherences
-subplot(4,4,2); imagesc(z_2, t_2, r12_re_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Real(rho_{12})')
-title('Real(rho_{12})')
-subplot(4,4,3); imagesc(z_2, t_2, r13_re_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Real(rho_{13})')
-title('Real(rho_{13})')
-subplot(4,4,4); imagesc(z_2, t_2, r14_re_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Real(rho_{14})')
-title('Real(rho_{14})')
-subplot(4,4,7); imagesc(z_2, t_2, r23_re_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Real(rho_{23})')
-title('Real(rho_{23})')
-subplot(4,4,8); imagesc(z_2, t_2, r24_re_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Real(rho_{24})')
-title('Real(rho_{24})')
-subplot(4,4,12); imagesc(z_2, t_2, r34_re_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Real(rho_{34})')
-title('Real(rho_{34})')
-% imaginary parts of coherences
-subplot(4,4,5); imagesc(z_2, t_2, r12_im_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Imag(rho_{12})')
-title('Imag(rho_{12})')
-subplot(4,4,9); imagesc(z_2, t_2, r13_im_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Imag(rho_{13})')
-title('Imag(rho_{13})')
-subplot(4,4,10); imagesc(z_2, t_2, r23_im_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Imag(rho_{23})')
-title('Imag(rho_{23})')
-subplot(4,4,13); imagesc(z_2, t_2, r14_im_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Imag(rho_{14})')
-title('Imag(rho_{14})')
-subplot(4,4,14); imagesc(z_2, t_2, r24_im_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Imag(rho_{24})')
-title('Imag(rho_{24})')
-subplot(4,4,15); imagesc(z_2, t_2, r34_im_out_2); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('Imag(rho_{34})')
-title('Imag(rho_{34})')
-
-
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp_I2.m
deleted file mode 100644
index 45d4913..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp_I2.m
+++ /dev/null
@@ -1,71 +0,0 @@
-Nlevels_no_dopler_with_z_4wm
-
-%% field propagation
-z_1=z_1*100; % z in cm
-t_1=t_1*1e6; % time now measured in uS
-figure(1)
-%set(gca,'fontsize',20); 
-imagesc(z_1, t_1, I2_out_1); colorbar
-tmin=-0.4;
-tmax= 0.4;
-ylim([tmin,tmax],'manual');	
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_2')
-title('I_2')
-
-xskip=1;
-yskip=10;
-%map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
-
-
-
-print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I2.eps')
-
-
-
-%% fields before and after the cell
-figure(2)
-%set(gca,'fontsize',30); 
-plot( ...
-	t_1,I2_out_1(:,1),'.-;before;', "linewidth", 4,  ...
-	t_1,I2_out_1(:,end), '-;after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I_2 (1/s)^2')
-title('I_2 before and after cell')
-legend('location', 'southwest');
-
-[b, a]=butter(3, 0.05);
-I2_out_after=I2_out_1(:,end);
-I2_out_after_filtered=filtfilt(b,a,I2_out_after);
-settling_time=0.8; %uS
-t_good_indx=t_1> min(t_1 + settling_time);
-[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
-delay_time=t_1(max_pos_after)-t_1(max_pos_before);
-printf('Second field delay time  = %f uS\n',delay_time);
-
-%set(gca,'fontsize',40); 
-%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell_I2.eps')
-
-figure(4)
-I2_max_in=max(I2_out_1(t_good_indx,1));
-I2_max_out=max(I2_out_1(t_good_indx,end));
-I2_in_norm=(I2_out_1(:,1))/I2_max_in;
-I2_out_norm=(I2_out_1(:,end))/I2_max_out;
-tmin=-0.05;
-tmax=0.05;
-indx=(t_1>=tmin & t_1<=tmax);  % soom in in time to this region
-plot( ...
-	t_1(indx),I2_in_norm(indx),'.-;before;', "linewidth", 4,  ...
-	t_1(indx),I2_out_norm(indx), '-;after;', "linewidth", 4 ...
-	)
-legend('location', 'southeast');
-xlim([tmin,tmax],'manual');	
-xlabel('t (uS)')
-ylabel('I_2')
-title('I_2 before and after cell normalized')
-%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','-depsc2', '-tight','-S200,120', 'probe_before_after_cell_I2_normalized.eps')
-
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/map2dat.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/map2dat.m
deleted file mode 120000
index 14fae30..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/map2dat.m
+++ /dev/null
@@ -1 +0,0 @@
-../map2dat.m
\ No newline at end of file
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/plot_fields_propagation_I2.gp b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/plot_fields_propagation_I2.gp
deleted file mode 100644
index a0cd31f..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/plot_fields_propagation_I2.gp
+++ /dev/null
@@ -1,15 +0,0 @@
-set terminal postscript portrait enhanced color solid size 5,3.5
-set output 'fields_propagation_I2.eps'
-set dgrid3d 100,100 qnorm 4
-set pm3d map
-#set contour
-set hidden3d
-set palette rgb 10,13,31 negative
-
-
-set xlabel "z (cm)"
-set ylabel "t ({/Symbol m}S)"
-set zlabel "I_2 (1/S)"
-set nokey
-#set view map
-splot [0:][-0.2:0.4] 'I2.dat' 
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/pp_I2.m
deleted file mode 100644
index ecfdbbb..0000000
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pulse_split/pp_I2.m
+++ /dev/null
@@ -1,68 +0,0 @@
-Nlevels_no_dopler_with_z_4wm
-
-%% field propagation
-z_1=z_1*100; % z in cm
-t_1=t_1*1e6; % time now measured in uS
-figure(1)
-set(gca,'fontsize',20); 
-imagesc(z_1, t_1, I2_out_1); colorbar
-xlabel('z (cm)')
-ylabel('t (uS)')
-zlabel('I_2')
-title('I_2')
-
-xskip=1;
-yskip=10;
-map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
-
-
-
-print('-color','fields_propagation_I2.eps')
-
-
-
-%% fields before and after the cell
-figure(2)
-%set(gca,'fontsize',30); 
-plot( ...
-	t_1,I2_out_1(:,1),'.-;before;', "linewidth", 4,  ...
-	t_1,I2_out_1(:,end), '-;after;', "linewidth", 4 ...
-	)
-xlabel('t (uS)')
-ylabel('I_2 (1/s)^2')
-title('I_2 before and after cell')
-legend('location', 'northeast');
-
-[b, a]=butter(3, 0.05);
-I2_out_after=I2_out_1(:,end);
-I2_out_after_filtered=filtfilt(b,a,I2_out_after);
-settling_time=0.01; %uS
-t_good_indx=t_1> min(t_1 + settling_time);
-[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
-delay_time=t_1(max_pos_after)-t_1(max_pos_before);
-printf('Second field delay time  = %f uS\n',delay_time);
-
-%set(gca,'fontsize',40); 
-set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','fields_before_after_cell_I2.eps')
-
-figure(4)
-I2_max_in=max(I2_out_1(t_good_indx,1));
-I2_max_out=max(I2_out_1(t_good_indx,end));
-I2_in_norm=(I2_out_1(:,1))/I2_max_in;
-I2_out_norm=(I2_out_1(:,end))/I2_max_out;
-tmin=-0.05;
-tmax=0.05;
-indx=(t_1>=tmin & t_1<=tmax);  % soom in in time to this region
-plot( ...
-	t_1(indx),I2_in_norm(indx),'.-;before;', "linewidth", 4,  ...
-	t_1(indx),I2_out_norm(indx), '-;after;', "linewidth", 4 ...
-	)
-xlim([tmin,tmax],'manual');	
-xlabel('t (uS)')
-ylabel('I_2')
-title('I_2 before and after cell normalized')
-legend('location', 'northeast');
-set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
-print('-color','probe_before_after_cell_I2_normalized.eps')
-
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_with_perturbations/Nlevels_no_dopler_with_z_4wm_with_perturbations.xmds b/xmds2/Nlevels_no_dopler_with_z_4wm_with_perturbations/Nlevels_no_dopler_with_z_4wm_with_perturbations.xmds
index e0440ef..e336db7 100644
--- a/xmds2/Nlevels_no_dopler_with_z_4wm_with_perturbations/Nlevels_no_dopler_with_z_4wm_with_perturbations.xmds
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_with_perturbations/Nlevels_no_dopler_with_z_4wm_with_perturbations.xmds
@@ -41,7 +41,7 @@
 
 		VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to
 		normalize it to something else look drho/dt equation.
-		No need to renormalizes eta as long as its express through i
+		No need to renormalizes eta as long as its express through
 		the upper level decay rate in the same units as Rabi frequency.
 	</description>
 
@@ -54,7 +54,8 @@
 				const double lambda=794.7e-9; //wavelength in m
 				const double Kvec = 2*M_PI/lambda;  // k-vector
 				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
-				const double eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				// eta will be calculated in the <arguments> section
+				double eta = 0;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
 				
 				//  ---------  Atom and cell properties -------------------------
 				// range of Maxwell distribution atomic velocities
@@ -62,9 +63,10 @@
 				// above mass expression is written as (expression is isotopic_mass * atomic_mass_unit)
 
 				// Average sqrt(v^2) in Maxwell distribution for one dimension
-				const double v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell related parameters will be calculated in <arguments> section
+				double v_thermal_averaged=0;
 				// Maxwell distribution velocities range to take in account in [m/s]
-				const double V_maxwell_min = -4*v_thermal_averaged, V_maxwell_max = -V_maxwell_min; // there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
+				double V_maxwell_min = 0, V_maxwell_max = 0;
 
 				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
 				const double gt=0.01 *(2*M_PI*1e6);
@@ -107,6 +109,18 @@
 					 I am guessing detunings are too large and thus it became a stiff equation-->
 			<!--! make sure it is not equal to zero!-->
 			<argument name="Temperature"  type="real" default_value="5" />
+			<!-- This will be executed after arguments/parameters are parsed -->
+			<!-- Read the code Luke: took me a while of reading the xmds2 sources to find it -->
+			<![CDATA[
+				// Average sqrt(v^2) in Maxwell distribution for one dimension
+				v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell distribution velocities range to take in account in [m/s]
+				// there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
+				V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; 
+
+				// eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;
+			]]>
     </arguments>
 		<bing />
 		<fftw plan="patient" />
diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds
index d63704f..bcb6b39 100644
--- a/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds
+++ b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds
@@ -37,7 +37,7 @@
 
 		VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to
 		normalize it to something else look drho/dt equation.
-		No need to renormalizes eta as long as its express through 
+		No need to renormalizes eta as long as its express through
 		the upper level decay rate in the same units as Rabi frequency.
 	</description>
 
@@ -50,7 +50,8 @@
 				const double lambda=794.7e-9; //wavelength in m
 				const double Kvec = 2*M_PI/lambda;  // k-vector
 				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
-				const double eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				// eta will be calculated in the <arguments> section
+				double eta = 0;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
 				
 				//  ---------  Atom and cell properties -------------------------
 				// range of Maxwell distribution atomic velocities
@@ -58,9 +59,10 @@
 				// above mass expression is written as (expression is isotopic_mass * atomic_mass_unit)
 
 				// Average sqrt(v^2) in Maxwell distribution for one dimension
-				const double v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell related parameters will be calculated in <arguments> section
+				double v_thermal_averaged=0;
 				// Maxwell distribution velocities range to take in account in [m/s]
-				const double V_maxwell_min = -4*v_thermal_averaged, V_maxwell_max = -V_maxwell_min; // there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
+				double V_maxwell_min = 0, V_maxwell_max = 0;
 
 				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
 				const double gt=0.01 *(2*M_PI*1e6);
@@ -106,6 +108,20 @@
 					 I am guessing detunings are too large and thus it became a stiff equation-->
 			<!--! make sure it is not equal to zero!-->
 			<argument name="Temperature"  type="real" default_value="5" />
+			<!-- This will be executed after arguments/parameters are parsed -->
+			<!-- Read the code Luke: took me a while of reading the xmds2 sources to find it -->
+			<![CDATA[
+				// Average sqrt(v^2) in Maxwell distribution for one dimension
+				if (Temperature == 0)
+					_LOG(_ERROR_LOG_LEVEL, "ERROR: Temperature should be >0 to provide range for Maxwell velocity distribution\n");
+				v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell distribution velocities range to take in account in [m/s]
+				// there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
+				V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; 
+
+				// eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;
+			]]>
     </arguments>
 		<bing />
 		<fftw plan="patient" />
@@ -117,7 +133,7 @@
 	<geometry>
 		<propagation_dimension> z </propagation_dimension>
 		<transverse_dimensions>
-			<dimension name="t"   lattice="1000"   domain="(-0.2e-6, 0.4e-6)" />
+			<dimension name="t"   lattice="1000"   domain="(-1e-6, 1e-6)" />
 			<dimension name="v"   lattice="100"   domain="(V_maxwell_min, V_maxwell_max)" />
 		</transverse_dimensions>
 	</geometry>
@@ -368,7 +384,7 @@
 	<!-- The output to generate -->
 	<output format="binary" filename="Nlevels_with_doppler_with_z_4wm.xsil">
 		<group>
-      <sampling basis="t(100) " initial_sample="yes">
+      <sampling basis="t(1000) " initial_sample="yes">
 				<dependencies>E_field_avgd</dependencies>
 				<moments>I1_out I2_out I3_out I4_out</moments>
 				<![CDATA[
diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/Makefile b/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/Makefile
index 11c0320..7dbc3f4 100644
--- a/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/Makefile
+++ b/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/Makefile
@@ -10,7 +10,15 @@ XSIL2GRAPHICS = xsil2graphics
 
 # fast light
 #PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=3.8e7 --E4o=1e1
-PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=0.1e7 --E2o=1e4 --E3o=.3e7 --E4o=0
+#PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=0.1e7 --E2o=1e4 --E3o=.3e7 --E4o=0 --Lcell=1.5e-2 --Temperature=.0001
+PARAMS =  \
+	  --Ndens=1e15 \
+	  --Lcell=10.0e-2 \
+	  --Temperature=5 \
+	  --Pwidth=0.4e-6 \
+	  --delta1=0 --delta2=0 --delta3=0 \
+	  --E1o=2e7 --E2o=1e2 --E3o=4e7 --E4o=1e0
+
 
 # slow light EIT
 #PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=0 --E4o=0
@@ -23,11 +31,13 @@ Nlevels_with_doppler_with_z_4wm.xsil: ../Nlevels_with_doppler_with_z_4wm.run
 %.m: %.xsil
 	$(XSIL2GRAPHICS) $<
 
-plot: $(M_FILES) $(GNUPLOT_FILES) 
-	octave pp_I2.m
+pretty_plots: plot $(GNUPLOT_FILES)
 	gnuplot plot_fields_propagation_I2.gp
 	gnuplot plot_fields_propagation_I4.gp
 
+plot: $(M_FILES)
+	octave pp_I2.m
+
 clean:
 	rm -f $(CC_FILES) $(RUN_FILES) $(M_FILES) $(XSIL_FILES) *.wisdom.fftw3 *.dat octave-core *.wisdom *.pdf
 	rm -f $(pdf_targets)
diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/README b/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/README
index 035178c..34ec5be 100644
--- a/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/README
+++ b/xmds2/Nlevels_with_doppler_with_z_4wm/fast_to_slow_switch/README
@@ -1,4 +1,13 @@
-Temperature = 1e-2K
+Fast light with advancement 0.006000 uS it is barely visible but there
+PARAMS =  \
+	  --Ndens=1e15 \
+	  --Lcell=10.0e-2 \
+	  --Temperature=5 \
+	  --Pwidth=0.4e-6 \
+	  --delta1=0 --delta2=0 --delta3=0 \
+	  --E1o=2e7 --E2o=1e2 --E3o=4e7 --E4o=1e0
+
+
 
 Fast light with decay in 1.5 cm cell
 PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=0.1e7 --E2o=1e4 --E3o=.3e7 --E4o=0
diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/pp_I2.m b/xmds2/Nlevels_with_doppler_with_z_4wm/pp_I2.m
index f632d0a..6b90568 100644
--- a/xmds2/Nlevels_with_doppler_with_z_4wm/pp_I2.m
+++ b/xmds2/Nlevels_with_doppler_with_z_4wm/pp_I2.m
@@ -31,8 +31,8 @@ tmax= 0.4;
 ylim([tmin,tmax],'manual');	
 xlabel('z (cm)')
 ylabel('t (uS)')
-zlabel('I_2')
-title('I_2')
+zlabel('I_4')
+title('I_4')
 print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I4.eps')
 map2dat('I4.dat',z_1,t_1, I4_out_1, xskip, yskip);
 
diff --git a/xmds2/TODO b/xmds2/TODO
new file mode 100644
index 0000000..c77e332
--- /dev/null
+++ b/xmds2/TODO
@@ -0,0 +1,3 @@
+- TODO Fix  Temperature > 10 misbehavior
+	Suspected stiff equation at this condition
+