delete some previous attempts of realistic Rb

1 files changed, 0 insertions, 642 deletions

diff --git a/xmds2/realistic_Rb/realistic_Rb.xmds b/xmds2/realistic_Rb/realistic_Rb.xmds
deleted file mode 100644
index 598715f..0000000
--- a/xmds2/realistic_Rb/realistic_Rb.xmds
+++ /dev/null
@@ -1,642 +0,0 @@
-<?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[
-				// Some numerical constants
-				const double pi = M_PI; 
-				// proportional to splitting ratios  sqrt(6) , sqrt(3), sqrt(2)
-				const double rt6 = 2.449489742783178;
-				const double rt3 = 1.7320508075688772;
-				const double rt2 = 1.4142135623730951;
-
-
-				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
-				// Fields k-vector
-				const double Kvec = 2*M_PI/lambda;  
-				// Simplified k-vectors
-				const double Kvec1 = Kvec, Kvec2=Kvec, Kvec3=Kvec;
-
-				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)]
-				double eta1=0, eta2=0, eta3=0;
-				
-				//  ---------  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 g1 = 3.612847284945266e7;
-				const double g2 = 3.8117309832741246e7;
-
-				// levels energy
-				const double ha0 = 2.1471788680034824e10;
-				const double ha1 = 2.558764384495815e9;
-				const double ha2 = 5.323020344462938e8;
-				const double hb2 = 7.85178251911697e7;
-
-				// Larmor frequency 
-				double WL=0;
-
-
-
-				complex  E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies
-
-				// density matrix elements which calculated via Hermitian property r_ij=conj(r_ji)
-				complex  
-						r1301, 
-						r1402, 
-						r0903, 
-						r1503, 
-						r1004, 
-						r1604, 
-						r1105, 
-						r0206, 
-						r1406, 
-						r0307, 
-						r0907, 
-						r1507, 
-						r0408, 
-						r1008, 
-						r1608, 
-						r1509, 
-						r1610; 
-
-
-				// 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;
-				// !FIXME over simplification: we should use relevant levels linewidths
-				eta1 = eta;
-				eta2 = eta;
-				eta3 = eta;
-			]]>
-    </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>
-				r0101a
-				r0113a
-				r0202a
-				r0214a
-				r0303a
-				r0309a
-				r0315a
-				r0404a
-				r0410a
-				r0416a
-				r0505a
-				r0511a
-				r0602a
-				r0606a
-				r0614a
-				r0703a
-				r0707a
-				r0709a
-				r0715a
-				r0804a
-				r0808a
-				r0810a
-				r0816a
-				r0909a
-				r0915a
-				r1010a
-				r1016a
-				r1111a
-				r1313a
-				r1414a
-				r1515a
-				r1616a
-		</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;
-
-				r0101a = r0101*prob_v_normalized;
-				r0113a = r0113*prob_v_normalized;
-				r0202a = r0202*prob_v_normalized;
-				r0214a = r0214*prob_v_normalized;
-				r0303a = r0303*prob_v_normalized;
-				r0309a = r0309*prob_v_normalized;
-				r0315a = r0315*prob_v_normalized;
-				r0404a = r0404*prob_v_normalized;
-				r0410a = r0410*prob_v_normalized;
-				r0416a = r0416*prob_v_normalized;
-				r0505a = r0505*prob_v_normalized;
-				r0511a = r0511*prob_v_normalized;
-				r0602a = r0602*prob_v_normalized;
-				r0606a = r0606*prob_v_normalized;
-				r0614a = r0614*prob_v_normalized;
-				r0703a = r0703*prob_v_normalized;
-				r0707a = r0707*prob_v_normalized;
-				r0709a = r0709*prob_v_normalized;
-				r0715a = r0715*prob_v_normalized;
-				r0804a = r0804*prob_v_normalized;
-				r0808a = r0808*prob_v_normalized;
-				r0810a = r0810*prob_v_normalized;
-				r0816a = r0816*prob_v_normalized;
-				r0909a = r0909*prob_v_normalized;
-				r0915a = r0915*prob_v_normalized;
-				r1010a = r1010*prob_v_normalized;
-				r1016a = r1016*prob_v_normalized;
-				r1111a = r1111*prob_v_normalized;
-				r1313a = r1313*prob_v_normalized;
-				r1414a = r1414*prob_v_normalized;
-				r1515a = r1515*prob_v_normalized;
-				r1616a = r1616*prob_v_normalized;
-			]]>
-		</evaluation>
-	</computed_vector>
-
-
-	<vector name="density_matrix" type="complex" initial_space="t">
-		<components>
-				r0101
-				r0113
-				r0202
-				r0214
-				r0303
-				r0309
-				r0315
-				r0404
-				r0410
-				r0416
-				r0505
-				r0511
-				r0602
-				r0606
-				r0614
-				r0703
-				r0707
-				r0709
-				r0715
-				r0804
-				r0808
-				r0810
-				r0816
-				r0909
-				r0915
-				r1010
-				r1016
-				r1111
-				r1313
-				r1414
-				r1515
-				r1616
-		</components>
-		<initialisation>
-			<!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
-			<![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
-				r0101 = 0.125;
-				r0113 = 0;
-				r0202 = 0.125;
-				r0214 = 0;
-				r0303 = 0.125;
-				r0309 = 0;
-				r0315 = 0;
-				r0404 = 0.125;
-				r0410 = 0;
-				r0416 = 0;
-				r0505 = 0.125;
-				r0511 = 0;
-				r0602 = 0;
-				r0606 = 0.125;
-				r0614 = 0;
-				r0703 = 0;
-				r0707 = 0.125;
-				r0709 = 0;
-				r0715 = 0;
-				r0804 = 0;
-				r0808 = 0.125;
-				r0810 = 0;
-				r0816 = 0;
-				r0909 = 0;
-				r0915 = 0;
-				r1010 = 0;
-				r1016 = 0;
-				r1111 = 0;
-				r1313 = 0;
-				r1414 = 0;
-				r1515 = 0;
-				r1616 = 0;
-			]]>
-		</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</samples>
-			<!--<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);
-
-						// Density matrix is Hermitian so we use r_ij=conj(r_ji)
-
-						r1301 = conj(r0113);
-						r1402 = conj(r0214);
-						r0903 = conj(r0309);
-						r1503 = conj(r0315);
-						r1004 = conj(r0410);
-						r1604 = conj(r0416);
-						r1105 = conj(r0511);
-						r0206 = conj(r0602);
-						r1406 = conj(r0614);
-						r0307 = conj(r0703);
-						r0907 = conj(r0709);
-						r1507 = conj(r0715);
-						r0408 = conj(r0804);
-						r1008 = conj(r0810);
-						r1608 = conj(r0816);
-						r1509 = conj(r0915);
-						r1610 = conj(r1016);
-
-						// Equations of motions according to Simon's mathematica code
-						dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1313)/6. - i*((r0113*E4a)/(4.*rt6) - (r1301*E4ac)/(4.*rt6));
-						dr0113_dt = (-(gt*r0113) - (gt + g2)*r0113)/2. - i*(WL*r0113 - ((2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0113 + (r0101*E4ac)/(4.*rt6) - (r1313*E4ac)/(4.*rt6));
-						dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (r1402*E4ac)/8.);
-						dr0214_dt = (-(gt*r0214) - (gt + g2)*r0214)/2. - i*((WL*r0214)/2. - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0214 - (r0206*E3ac)/(8.*rt3) + (r0202*E4ac)/8. - (r1414*E4ac)/8.);
-						dr0303_dt = gt/8. - gt*r0303 + (g1*r0909)/12. + (g1*r1010)/3. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1515)/4. - i*((r0309*E1a)/(4.*rt6) + (r0315*E4a)/8. - (r0903*E1ac)/(4.*rt6) - (r1503*E4ac)/8.);
-						dr0309_dt = (-(gt*r0309) - (gt + g1)*r0309)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0309) + (r0303*E1ac)/(4.*rt6) - (r0909*E1ac)/(4.*rt6) - (r0307*E2ac)/(4.*rt6) - (r1509*E4ac)/8.);
-						dr0315_dt = (-(gt*r0315) - (gt + g2)*r0315)/2. - i*(-(((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0315) - (r0915*E1ac)/(4.*rt6) - (r0307*E3ac)/8. + (r0303*E4ac)/8. - (r1515*E4ac)/8.);
-						dr0404_dt = gt/8. - gt*r0404 + (g1*r1010)/4. + (g1*r1111)/4. + (g2*r1414)/4. + (g2*r1515)/12. + (g2*r1616)/6. - i*((r0410*E1a)/(4.*rt2) + (r0416*E4a)/(4.*rt6) - (r1004*E1ac)/(4.*rt2) - (r1604*E4ac)/(4.*rt6));
-						dr0410_dt = (-(gt*r0410) - (gt + g1)*r0410)/2. - i*(-(WL*r0410)/2. + (delta1 + v*Kvec1)*r0410 + (r0404*E1ac)/(4.*rt2) - (r1010*E1ac)/(4.*rt2) - (r0408*E2ac)/(4.*rt6) - (r1610*E4ac)/(4.*rt6));
-						dr0416_dt = (-(gt*r0416) - (gt + g2)*r0416)/2. - i*(-(WL*r0416)/2. - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0416 - (r1016*E1ac)/(4.*rt2) - (r0408*E3ac)/(4.*rt2) + (r0404*E4ac)/(4.*rt6) - (r1616*E4ac)/(4.*rt6));
-						dr0505_dt = gt/8. - gt*r0505 + (g1*r1111)/2. + (g2*r1515)/6. + (g2*r1616)/3. - i*((r0511*E1a)/4. - (r1105*E1ac)/4.);
-						dr0511_dt = (-(gt*r0511) - (gt + g1)*r0511)/2. - i*(-(WL*r0511) - (WL/6. - delta1 - v*Kvec1)*r0511 + (r0505*E1ac)/4. - (r1111*E1ac)/4.);
-						dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (r1402*E3ac)/(8.*rt3));
-						dr0606_dt = gt/8. - gt*r0606 + (g1*r0909)/12. + (g1*r1010)/12. + (g2*r1313)/4. + (g2*r1414)/12. - i*(-(r0614*E3a)/(8.*rt3) + (r1406*E3ac)/(8.*rt3));
-						dr0614_dt = (-(gt*r0614) - (gt + g2)*r0614)/2. - i*((-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0614 - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0614 - (r0606*E3ac)/(8.*rt3) + (r1414*E3ac)/(8.*rt3) + (r0602*E4ac)/8.);
-						dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (r0903*E2ac)/(4.*rt6) + (r1503*E3ac)/8.);
-						dr0707_dt = gt/8. - gt*r0707 + (g1*r0909)/12. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1414)/3. + (g2*r1515)/4. - i*(-(r0709*E2a)/(4.*rt6) - (r0715*E3a)/8. + (r0907*E2ac)/(4.*rt6) + (r1507*E3ac)/8.);
-						dr0709_dt = (-(gt*r0709) - (gt + g1)*r0709)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0709) + (-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0709 + (r0703*E1ac)/(4.*rt6) - (r0707*E2ac)/(4.*rt6) + (r0909*E2ac)/(4.*rt6) + (r1509*E3ac)/8.);
-						dr0715_dt = (-(gt*r0715) - (gt + g2)*r0715)/2. - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0715 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0715 + (r0915*E2ac)/(4.*rt6) - (r0707*E3ac)/8. + (r1515*E3ac)/8. + (r0703*E4ac)/8.);
-						dr0804_dt = -(gt*r0804) - i*((WL*r0804)/2. + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0804 + (r0810*E1a)/(4.*rt2) + (r0816*E4a)/(4.*rt6) + (r1004*E2ac)/(4.*rt6) + (r1604*E3ac)/(4.*rt2));
-						dr0808_dt = gt/8. - gt*r0808 + (g1*r1010)/12. + (g1*r1111)/12. + (g2*r1414)/12. + (g2*r1515)/4. + (g2*r1616)/2. - i*(-(r0810*E2a)/(4.*rt6) - (r0816*E3a)/(4.*rt2) + (r1008*E2ac)/(4.*rt6) + (r1608*E3ac)/(4.*rt2));
-						dr0810_dt = (-(gt*r0810) - (gt + g1)*r0810)/2. - i*((delta1 + v*Kvec1)*r0810 + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0810 + (r0804*E1ac)/(4.*rt2) - (r0808*E2ac)/(4.*rt6) + (r1010*E2ac)/(4.*rt6) + (r1610*E3ac)/(4.*rt2));
-						dr0816_dt = (-(gt*r0816) - (gt + g2)*r0816)/2. - i*((WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0816 - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0816 + (r1016*E2ac)/(4.*rt6) - (r0808*E3ac)/(4.*rt2) + (r1616*E3ac)/(4.*rt2) + (r0804*E4ac)/(4.*rt6));
-						dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (r0903*E1ac)/(4.*rt6) - (r0907*E2ac)/(4.*rt6));
-						dr0915_dt = (-((gt + g1)*r0915) - (gt + g2)*r0915)/2. - i*((-WL/6. - delta1 - v*Kvec1)*r0915 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0915 - (r0315*E1a)/(4.*rt6) + (r0715*E2a)/(4.*rt6) - (r0907*E3ac)/8. + (r0903*E4ac)/8.);
-						dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (r1004*E1ac)/(4.*rt2) - (r1008*E2ac)/(4.*rt6));
-						dr1016_dt = (-((gt + g1)*r1016) - (gt + g2)*r1016)/2. - i*(-((delta1 + v*Kvec1)*r1016) - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r1016 - (r0416*E1a)/(4.*rt2) + (r0816*E2a)/(4.*rt6) - (r1008*E3ac)/(4.*rt2) + (r1004*E4ac)/(4.*rt6));
-						dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (r1105*E1ac)/4.);
-						dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (r1301*E4ac)/(4.*rt6));
-						dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (r1406*E3ac)/(8.*rt3) + (r1402*E4ac)/8.);
-						dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (r1507*E3ac)/8. + (r1503*E4ac)/8.);
-						dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (r1608*E3ac)/(4.*rt2) + (r1604*E4ac)/(4.*rt6));
-					]]>
-        </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 = 0.16666666666666666*eta1*(2.449489742783178*r0309 + 4.242640687119286*r0410 + 6.*r0511) - Lt[E1];
-					dE2_dz = -0.8164965809277261*eta1*(r0709 + r0810) - Lt[E2];
-					dE3_dz = -1.*eta2*(1.7320508075688772*r0614 + 3.*r0715 + 4.242640687119286*r0816) - Lt[E3];
-					dE4_dz = (4*eta2*(2.449489742783178*r0113 + 3*r0214 + 3*r0315 + 2.449489742783178*r0416))/3. - 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: 
--->