changed some substitution to better match our XMDS

5 files changed, 699 insertions, 0 deletions

diff --git a/xmds2/realistic_Rb/Makefile b/xmds2/realistic_Rb/Makefile
new file mode 100644
index 0000000..e088e15
--- /dev/null
+++ b/xmds2/realistic_Rb/Makefile
@@ -0,0 +1,55 @@
+###  -*- 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: $(RUN_FILES) 
+
+%.run: %.xmds
+	$(XMDS) $<
+	mv $(patsubst %.xmds,%,$<) $@
+
+%.xsil: %.run
+	./$<  --E1o=0 --E3o=0
+
+%.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)
+	rm -f $(eps_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/realistic_Rb/Makefile.fig b/xmds2/realistic_Rb/Makefile.fig
new file mode 100644
index 0000000..71fba31
--- /dev/null
+++ b/xmds2/realistic_Rb/Makefile.fig
@@ -0,0 +1,39 @@
+###  -*- make -*-
+### This makefile can be used to build and run the XMDS examples
+
+GNUPLOT_FILES = $(wildcard *.gp)
+
+SCRIPTS_DIR = .
+
+
+eps_targets  = $(wildcard *.eps)
+eps_bb_targets  =  $(eps_targets:%.eps=%.bb)
+pdf_targets  = $(eps_targets:%.eps=%.pdf)	
+png_targets  = $(pdf_targets:%.pdf=%.png)
+
+fig: png
+
+pdf: $(pdf_targets) 
+
+# fix bounding box
+$(eps_bb_targets): %.bb : %.eps
+	@cat $< | ps2eps -q -B > $@
+
+$(pdf_targets): %.pdf : %.bb
+	epspdf $< $@
+
+png: $(png_targets) 
+
+$(png_targets): %.png : %.pdf
+	convert -density 300 $< $@
+
+clean:
+	rm -f *.pdf
+	rm -f *.eps
+	rm -f *.bb
+
+real_clean: clean
+	rm -f *.png
+
+.PRECIOUS: %.run %.xsil %.m
+.PHONY: all clean
diff --git a/xmds2/realistic_Rb/Makefile.par b/xmds2/realistic_Rb/Makefile.par
new file mode 100644
index 0000000..cc21c7b
--- /dev/null
+++ b/xmds2/realistic_Rb/Makefile.par
@@ -0,0 +1,32 @@
+###  -*- 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
+
+
+PARAMS_FILES = $(wildcard *.params)
+PP_DIR = $(PARAMS_FILES:%.params=%)
+CALC_XSIL_FILES = $(PARAMS_FILES:%.params=%/data.xsil)
+
+default: $(CALC_XSIL_FILES) 
+
+$(CALC_XSIL_FILES): %/data.xsil : % %.params  
+	echo processing $$(dirname $(@)) dir
+	$(MAKE) -C $$(dirname $(@)) -f ../Makefile.pp SCRIPTS_DIR=../ PARAMS_FILE=../$<.params  
+
+$(PP_DIR): % : %.params
+	echo need to make dir
+	[ -d $@ ] ||  mkdir -p  $@
+
+clean:
+	for d in $(PP_DIR); \
+	do $(MAKE) -C $$d SCRIPTS_DIR=../ -f ../Makefile.pp $@; done
+
+real_clean:
+	for d in $(PP_DIR); \
+	do $(MAKE) -C $$d SCRIPTS_DIR=../ -f ../Makefile.pp $@; done
+
+.PHONY: all clean
diff --git a/xmds2/realistic_Rb/Makefile.pp b/xmds2/realistic_Rb/Makefile.pp
new file mode 100644
index 0000000..83e5bec
--- /dev/null
+++ b/xmds2/realistic_Rb/Makefile.pp
@@ -0,0 +1,62 @@
+###  -*- make -*-
+### This makefile can be used to build and run the XMDS examples
+
+
+XSIL_FILES = Nlevels_with_doppler_with_z_4wm.xsil
+M_FILES = $(patsubst %.xsil,%.m,$(XSIL_FILES))
+GNUPLOT_FILES = $(wildcard *.gp)
+
+SCRIPTS_DIR = .
+
+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 --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
+
+
+include $(PARAMS_FILE)
+
+
+# slow light EIT
+#PARAMS =  --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=0 --E4o=0
+
+all: $(XSIL_FILES) Nlevels_with_doppler_with_z_4wm.xsil $(M_FILES) plot fig
+
+Nlevels_with_doppler_with_z_4wm.xsil: ../Nlevels_with_doppler_with_z_4wm.run $(PARAMS_FILE)
+	$< $(PARAMS) | grep "Time elapsed for simulation is:" > exact_analysis_execution_time.txt
+
+%.m: %.xsil
+	$(XSIL2GRAPHICS) $<
+
+fig: pp_I2.stamp
+	$(MAKE) -f $(SCRIPTS_DIR)/Makefile.fig $@
+
+
+pretty_plots: pp_I2.stamp $(GNUPLOT_FILES)
+	gnuplot $(SCRIPTS_DIR)/plot_fields_propagation_I2.gp
+	gnuplot $(SCRIPTS_DIR)/plot_fields_propagation_I4.gp
+
+plot: pp_I2.stamp
+
+pp_I2.stamp: $(XSIL_FILES) $(M_FILES)
+	octave -q $(SCRIPTS_DIR)/pp_I2.m
+
+clean:
+	rm -f $(M_FILES) $(XSIL_FILES) *.dat octave-core 
+	rm -f pp_I2.stamp
+	$(MAKE) -f $(SCRIPTS_DIR)/Makefile.fig $@
+
+real_clean: clean
+	$(MAKE) -f $(SCRIPTS_DIR)/Makefile.fig $@
+	rm -f exact_analysis_execution_time.txt
+
+.PRECIOUS: %.run %.xsil %.m
+.PHONY: all clean
diff --git a/xmds2/realistic_Rb/realistic_Rb.xmds b/xmds2/realistic_Rb/realistic_Rb.xmds
new file mode 100644
index 0000000..6a40191
--- /dev/null
+++ b/xmds2/realistic_Rb/realistic_Rb.xmds
@@ -0,0 +1,511 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+	<name>realistic_Rb</name>
+
+	<author>Eugeniy Mikhailov</author>
+	<description>
+		License GPL.
+
+		Solving simplified Rb atom model
+		with fields propagation along spatial axis Z
+		with 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.
+
+
+		*              --------------- | F=1, 2P_3/2 >
+		*                \        \    
+		*                 \ E3_r   \    -------- | F=2, 2P_+1/2 >
+		*                  \   E4_r \   /    \
+		*                   \        \ / E2_l \ 
+		*                    \        /        \ E1_l
+		*   | F=2, 2S_1/2 > --------------      \
+		*                               \        \
+		*                                \        \
+		*                               ------------- | F=1, 2S_1/2 >
+		*            
+
+
+		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
+		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]
+				// 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
+				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
+				// 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]
+				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);
+				// 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  E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies
+
+				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
+
+				// inner use variables 
+				double probability_v; // will be used as p(v) in Maxwell distribution
+
+			]]>
+		</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" />
+			<!-- 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 />
+		<diagnostics /> 
+		<fftw plan="patient" threads="1" />
+		<!-- I don't see any speed up on 6 core CPU even if use threads="6" -->
+		<openmp />
+		<auto_vectorise />
+		<halt_non_finite />
+	</features>
+
+	<!-- 'z', 't', and 'v'  to have dimensions [m], [s], and [m/s]   -->
+	<geometry>
+		<propagation_dimension> z </propagation_dimension>
+		<transverse_dimensions>
+			<!-- IMPORTANT: looks like having a lot of points in time helps with convergence.
+					 I suspect that time spacing should be small enough to catch
+					 all pulse harmonics and more importantly 1/dt should be larger than
+					 the largest detuning (including Doppler shifts).
+					 Unfortunately calculation time is proportional to lattice size
+					 so we cannot just blindly increase it.
+					 Some rules of thumb:  
+						* lattice="1000"   domain="(-1e-6, 1e-6)" 
+							was good enough detunings up to 155 MHz (980 rad/s) notice that 1/dt=500 MHz
+						* lattice="10000"   domain="(-1e-6, 1e-6)" 
+							works for Doppler averaging in up to 400K for Rb when lasers are zero detuned
+			 -->
+			<dimension name="t"   lattice="10000"   domain="(-1e-6, 1e-6)" />
+			<dimension name="v"   lattice="100"   domain="(V_maxwell_min, V_maxwell_max)" />
+		</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>
+
+	<!--Maxwell distribution probability p(v)-->
+	<computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real">
+		<components>probability_v</components>
+		<evaluation>
+			<![CDATA[
+			// TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during  the script lifetime since 'v' is fixed.
+			probability_v=1.0/(v_thermal_averaged*sqrt(2*M_PI)) * exp( - mod2(v/v_thermal_averaged)/2.0 ); 
+			]]>
+		</evaluation>
+	</computed_vector>
+
+	<!--Maxwell distribution norm sum(p(v))
+			 Needed since we sum over the grid instead of true integral,
+			 we also have finite cut off velocities-->
+	<computed_vector name="Maxwell_distribution_probabilities_norm" dimensions="" type="real">
+		<components>probability_v_norm</components>
+		<evaluation>
+			<dependencies basis="v">Maxwell_distribution_probabilities</dependencies>
+			<![CDATA[
+			// TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during  the script lifetime since 'v' is fixed.
+			probability_v_norm=probability_v;
+			]]>
+		</evaluation>
+	</computed_vector>
+
+
+	<!-- Averaged across Maxwell distribution fields amplitudes -->
+	<computed_vector name="E_field_avgd" dimensions="t" type="complex">
+		<components>E1a E2a E3a E4a</components>
+		<evaluation>
+			<dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
+			<![CDATA[
+			double prob_v_normalized=probability_v/probability_v_norm;
+			E1a=E1*prob_v_normalized;
+			E2a=E2*prob_v_normalized;
+			E3a=E3*prob_v_normalized;
+			E4a=E4*prob_v_normalized;
+			]]>
+		</evaluation>
+	</computed_vector>
+
+	<!-- Averaged across Maxwell distribution density matrix components -->
+	<computed_vector name="density_matrix_averaged" dimensions="t" type="complex">
+		<components>r11a r22a r33a r12a r13a r14a r23a r24a r34a r44a</components>
+		<evaluation>
+			<dependencies basis="v">density_matrix Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
+			<![CDATA[
+			double prob_v_normalized=probability_v/probability_v_norm;
+			r11a=r11*prob_v_normalized;
+			r22a=r22*prob_v_normalized;
+			r33a=r33*prob_v_normalized;
+			r12a=r12*prob_v_normalized;
+			r13a=r13*prob_v_normalized;
+			r14a=r14*prob_v_normalized;
+			r23a=r23*prob_v_normalized;
+			r24a=r24*prob_v_normalized;
+			r34a=r34*prob_v_normalized;
+			r44a=r44*prob_v_normalized;
+			]]>
+		</evaluation>
+	</computed_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_avgd</dependencies>
+				<![CDATA[
+						E1ac = conj(E1a);
+						E2ac = conj(E2a);
+						E3ac = conj(E3a);
+						E4ac = conj(E4a);
+
+				// IMPORTANT: assumes no detunings
+				r11 = (gt*(4*mod2(E1a) + (G3 + gt)*(G3 + 2*gt))*((G4 +
+				2*gt)*(4*mod2(E3a) + gt*(G4 + gt)) + 4*mod2(E3a)*G4*(R41 -
+				R42)))/(2.*(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
+				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
+				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
+				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
+
+				r13 = -((E1ac*gt*(G3 + gt)*i*((G4 + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))
+				+ 4*mod2(E3a)*G4*(R41 -
+				R42)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
+				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
+				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
+				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
+
+				r22 = (gt*(4*mod2(E3a) + (G4 + gt)*(G4 + 2*gt))*((G3 +
+				2*gt)*(4*mod2(E1a) + gt*(G3 + gt)) + 4*mod2(E1a)*G3*(-R31 +
+				R32)))/(2.*(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
+				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
+				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
+				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
+
+				r24 = -((E3ac*gt*(G4 + gt)*i*((G3 + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))
+				+ 4*mod2(E1a)*G3*(-R31 +
+				R32)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
+				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
+				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
+				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
+
+				r33 = (2*mod2(E1a)*gt*((G4 + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt)) +
+				4*mod2(E3a)*G4*(R41 -
+				R42)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
+				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
+				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
+				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42));
+
+				r44 = (2*mod2(E3a)*gt*((G3 + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt)) +
+				4*mod2(E1a)*G3*(-R31 +
+				R32)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
+				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
+				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
+				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*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"> 
+		<!--<integrate algorithm="SI" steps="200" 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>100 100</samples>
+			<!--Use the next line for debuging to see velocity dependence. Uncomment/switch on output groups 3,4-->
+			<!--<samples>100 100 100 100</samples>--> 
+			<operators>
+        <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
+					<integration_vectors>density_matrix</integration_vectors>
+          <dependencies>E_field_avgd</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[
+						E1ac = conj(E1a);
+						E2ac = conj(E2a);
+						E3ac = conj(E3a);
+						E4ac = conj(E4a);
+
+						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*(-(E1a*r13) - E4a*r14 + E1ac*r31 + E4ac*r41) + G3*r33*R31 + G4*r44*R41;
+						dr12_dt = -(gt*r12) + i*((-(delta1 + Kvec*v) + (delta2 + Kvec*v))*r12 - E2a*r13 - E3a*r14 + E1ac*r32 + E4ac*r42);
+						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1ac*r11) - E2ac*r12 - (delta1 + Kvec*v)*r13 + E1ac*r33 + E4ac*r43);
+						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4ac*r11) - E3ac*r12 - ((delta1 + Kvec*v) - (delta2 + Kvec*v) + (delta3 + Kvec*v))*r14 + E1ac*r34 + E4ac*r44);
+						dr22_dt = gt/2 - gt*r22 + i*(-(E2a*r23) - E3a*r24 + E2ac*r32 + E3ac*r42) + G3*r33*R32 + G4*r44*R42;
+						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1ac*r21) - E2ac*r22 - (delta2 + Kvec*v)*r23 + E2ac*r33 + E3ac*r43);
+						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4ac*r21) - E3ac*r22 - (delta3 + Kvec*v)*r24 + E2ac*r34 + E3ac*r44);
+						dr33_dt = i*(E1a*r13 + E2a*r23 - E1ac*r31 - E2ac*r32) - (G3 + gt)*r33;
+
+						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1a*r14 + E2a*r24 - E4ac*r31 - E3ac*r32 + ((delta2 + Kvec*v) - (delta3 + Kvec*v))*r34);
+						dr44_dt = i*(E4a*r14 + E3a*r24 - E4ac*r41 - E3ac*r42) - (G4 + gt)*r44;
+					]]>
+        </operator>
+				<!--
+							 According to xmds2 docs operator kind="ip" should be faster
+							 but our codes runs about 5% to 10% slower with it.
+							 Maybe because we very close to the stiff condition so I use "ex" kind
+							 <operator kind="ip" constant="yes">
+					 -->
+				<operator kind="ex" constant="yes" type="imaginary">
+					<operator_names>Lt</operator_names>
+					<![CDATA[
+					Lt = -i/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="realistic_Rb.xsil">
+		<group>
+      <sampling basis="t(1000) " initial_sample="yes">
+				<dependencies>E_field_avgd</dependencies>
+				<moments>I1_out I2_out I3_out I4_out</moments>
+				<![CDATA[
+				I1_out = mod2(E1a);
+				I2_out = mod2(E2a);
+				I3_out = mod2(E3a);
+				I4_out = mod2(E4a);
+				]]>
+			</sampling>
+		</group>
+
+		<group>
+      <sampling basis="t(100) v(10)" initial_sample="yes">
+				<dependencies>density_matrix_averaged</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 = r11a.Re();
+				r22_out = r22a.Re();
+				r33_out = r33a.Re();
+				r44_out = r44a.Re();
+				// coherences output 
+				r12_re_out = r12a.Re();
+				r12_im_out = r12a.Im();
+				r13_re_out = r13a.Re();
+				r13_im_out = r13a.Im();
+				r14_re_out = r14a.Re();
+				r14_im_out = r14a.Im();
+				r23_re_out = r23a.Re();
+				r23_im_out = r23a.Im();
+				r24_re_out = r24a.Re();
+				r24_im_out = r24a.Im();
+				r34_re_out = r34a.Re();
+				r34_im_out = r34a.Im();
+				]]>
+			</sampling>
+		</group>
+
+		<!-- use the following two groups only for debuging 
+				 otherwise they are quite useless and have to much information 
+				 in 3D space (z,t,v) -->
+		<!--
+		<group>
+      <sampling basis="t(100) v(10)" initial_sample="yes">
+				<dependencies>E_field</dependencies>
+				<moments>I1_out_v I2_out_v I3_out_v I4_out_v</moments>
+				<![CDATA[
+				// light field intensity distribution in velocity subgroups 
+				I1_out_v = mod2(E1);
+				I2_out_v = mod2(E2);
+				I3_out_v = mod2(E3);
+				I4_out_v = mod2(E4);
+				]]>
+			</sampling>
+		</group>
+
+		<group>
+      <sampling basis="t(100) v(10)" initial_sample="yes">
+				<dependencies>density_matrix</dependencies>
+				<moments>
+					r11_out_v r22_out_v r33_out_v r44_out_v 
+					r12_re_out_v r12_im_out_v r13_re_out_v r13_im_out_v r14_re_out_v r14_im_out_v
+					                      r23_re_out_v r23_im_out_v r24_re_out_v r24_im_out_v 
+					                                            r34_re_out_v r34_im_out_v
+				</moments>
+				<![CDATA[
+				// density matrix distribution in velocity subgroups 
+				// populations output 
+				r11_out_v = r11.Re();
+				r22_out_v = r22.Re();
+				r33_out_v = r33.Re();
+				r44_out_v = r44.Re();
+				// coherences output 
+				r12_re_out_v = r12.Re();
+				r12_im_out_v = r12.Im();
+				r13_re_out_v = r13.Re();
+				r13_im_out_v = r13.Im();
+				r14_re_out_v = r14.Re();
+				r14_im_out_v = r14.Im();
+				r23_re_out_v = r23.Re();
+				r23_im_out_v = r23.Im();
+				r24_re_out_v = r24.Re();
+				r24_im_out_v = r24.Im();
+				r34_re_out_v = r34.Re();
+				r34_im_out_v = r34.Im();
+				]]>
+			</sampling>
+		</group>
+		-->
+
+	</output>
+
+</simulation>
+	
+<!-- 
+vim: ts=2 sw=2 foldmethod=indent: 
+-->