added test files for simulations

5 files changed, 1719 insertions, 0 deletions

diff --git a/xmds2/realistic_Rb/Makefile b/xmds2/realistic_Rb/Makefile
index e088e15..34d73e1 100644
--- a/xmds2/realistic_Rb/Makefile
+++ b/xmds2/realistic_Rb/Makefile
@@ -46,6 +46,9 @@ $(png_targets): %.png : %.pdf
 
 pdf: $(pdf_targets) 
 
+test:
+	cd tests; ./run_tests.py
+
 $(pdf_targets): %.pdf : %.eps
 	cat $< | ps2eps -B > __tt.eps
 	epspdf __tt.eps $@
diff --git a/xmds2/realistic_Rb/tests/run_tests.py b/xmds2/realistic_Rb/tests/run_tests.py
new file mode 100755
index 0000000..61f4eea
--- /dev/null
+++ b/xmds2/realistic_Rb/tests/run_tests.py
@@ -0,0 +1,350 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+run_tests.py
+
+Created by Graham Dennis on 2008-06-15.
+
+Copyright (c) 2008-2012, Graham Dennis
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+import xpdeint.Python24Support
+
+import os
+import re
+import sys
+import getopt
+import shutil
+import hashlib
+import unittest
+import subprocess
+
+from xml.dom import minidom
+import xpdeint.minidom_extras
+from xpdeint import CodeParser
+
+from xpdeint.XSILFile import XSILFile
+
+import numpy
+
+help_message = '''
+The help message goes here.
+'''
+
+
+class Usage(Exception):
+  def __init__(self, msg):
+    self.msg = msg
+
+def pass_nan_test(array1, array2):
+    """Return `True` if isNaN(`array1`) == isNaN(`array2`)"""
+    # NaN test. array2 is allowed to be NaN at an index if array1 is also NaN there.
+    nanTestPassed = numpy.equal(numpy.isnan(array1), numpy.isnan(array2)).all()
+    return nanTestPassed
+
+def array_approx_equal(array1, array2, absTol, relTol):
+  """Return `True` if all of (`array1` - `array2`) <= `absTol` or (`array1` - `array2`) <= `relTol` * `array2`"""
+  diff = array1-array2
+  # NaN values would fail this test. So we have to exclude them. But only exclude them if array2 (the expected results)
+  # have a NaN
+  return numpy.logical_or(numpy.logical_or(numpy.abs(diff) <= 0.5 * relTol * (numpy.abs(array2) + numpy.abs(array1)), numpy.abs(diff) <= absTol), numpy.isnan(array2)).all()
+
+def scriptTestingFunction(root, scriptName, testDir, absPath, self):
+  if not os.path.exists(testDir):
+    os.makedirs(testDir)
+  
+  proc = subprocess.Popen('xmds2 --no-version ' + absPath,
+                          shell=True,
+                          stdout=subprocess.PIPE,
+                          stderr=subprocess.PIPE,
+                          cwd=testDir)
+  (stdout, stderr) = proc.communicate()
+  returnCode = proc.wait()
+  
+  message = ''.join(["\n%(handleName)s:\n%(content)s" % locals() for handleName, content in [('stdout', stdout), ('stderr', stderr)] if content])
+  
+  # A few tests require XMDS1. If XMDS1 isn't present we should just
+  # skip that test rather than failing.
+  # The skip functionality for the unittest class is only available
+  # in python 2.7 and later, so check for that too.
+  if returnCode != 0 and sys.version_info[:2] >= (2, 7):
+    if re.search(r'^The missing \w+ feature\(s\) were: .*xmds.*', message, re.MULTILINE):
+      self.skipTest("Skipping test as XMDS1 is required and not installed")
+
+  # A few tests require specific features.  If it isn't available, skip the test
+  # rather than failing.
+  # The skip functionality for the unittest class is only available
+  # in python 2.7 and later, so check for that too.
+  if returnCode != 0 and sys.version_info[:2] >= (2, 7):
+    if re.search(r'^The missing \w+ feature\(s\) were:', message, re.MULTILINE):
+      self.skipTest("Skipping test as feature required is not installed")
+
+  self.assert_(returnCode == 0, ("Failed to compile." % locals()) + message)
+  
+  xmlDocument = minidom.parse(absPath)
+  simulationElement = xmlDocument.getChildElementByTagName('simulation')
+  nameElement = simulationElement.getChildElementByTagName('name')
+  testingElement = simulationElement.getChildElementByTagName('testing')
+  
+  simulationName = nameElement.innerText()
+  
+  # If the source is the same as the last known good, then we don't need to compile or execute the simulation.
+  sourceFilePath = os.path.join(testDir, simulationName + '.cc')
+  checksumFilePath = os.path.join(testDir, simulationName + '_last_known_good.checksum')
+  sourceContents = file(sourceFilePath).read()
+  h = hashlib.sha1()
+  h.update(sourceContents)
+  currentChecksum = h.hexdigest()
+  
+  if os.path.exists(checksumFilePath):
+    lastKnownGoodChecksum = file(checksumFilePath).read()
+    
+    if lastKnownGoodChecksum == currentChecksum:
+      # The checksums check out, so we don't need to go any further
+      return
+  
+  # Now we have compiled, we need to copy any input data needed and then run the simulation
+  inputXSILElements = testingElement.getChildElementsByTagName('input_xsil_file', optional=True)
+  
+  filesToCopy = []
+  
+  for inputXSILElement in inputXSILElements:
+    name = inputXSILElement.getAttribute('name').strip()
+    filesToCopy.append(name)
+    inputXSILFile = XSILFile(os.path.join(os.path.split(absPath)[0], name), loadData=False)
+    filesToCopy.extend([os.path.join(os.path.split(name)[0], xsil.data.filename) for xsil in inputXSILFile.xsilObjects if hasattr(xsil.data, 'filename')])
+  
+  for fileToCopy in filesToCopy:
+    sourceFile = os.path.join(os.path.split(absPath)[0], fileToCopy)
+    shutil.copy(sourceFile, testDir)
+  
+  # Allow command-line arguments to be specified for the simulation
+  commandLineElement = testingElement.getChildElementByTagName('command_line', optional=True)
+  argumentsElement = testingElement.getChildElementByTagName('arguments', optional=True)
+  commandLineString = './' + simulationName
+  if commandLineElement:
+    # The command line element overrides the prefix
+    commandLineString = commandLineElement.innerText().strip()
+  if argumentsElement:
+    commandLineString += ' ' + argumentsElement.innerText().strip()
+  
+  simulationProc = subprocess.Popen(commandLineString,
+                                    shell=True,
+                                    stdout=subprocess.PIPE,
+                                    stderr=subprocess.PIPE,
+                                    cwd=testDir)
+  (stdout, stderr) = simulationProc.communicate()
+  returnCode = simulationProc.wait()
+  
+  self.assert_(returnCode == 0, "Failed to execute compiled simulation correctly." % locals())
+  
+  # The next thing to check is that the generated data agrees with the expected data to within the set error margins.
+  xsilFileElements = testingElement.getChildElementsByTagName('xsil_file', optional=True)
+  for xsilFileElement in xsilFileElements:
+    sourceFile = xsilFileElement.getAttribute('name').strip()
+    expectedResultsFile = xsilFileElement.getAttribute('expected').strip()
+    # Defaults
+    absoluteTolerance = 0 
+    relativeTolerance = 1e-9
+    
+    if xsilFileElement.hasAttribute('absolute_tolerance'):
+      absoluteTolerance = float(xsilFileElement.getAttribute('absolute_tolerance'))
+    if xsilFileElement.hasAttribute('relative_tolerance'):
+      relativeTolerance = float(xsilFileElement.getAttribute('relative_tolerance'))
+    
+    resultsFullPath = os.path.join(testDir, sourceFile)
+    results = XSILFile(resultsFullPath)
+    expectedResultsFullPath = os.path.join(os.path.split(absPath)[0], expectedResultsFile)
+    if not os.path.exists(expectedResultsFullPath):
+      print >> sys.stderr, "Expected results file '%(expectedResultsFile)s' missing. Using current. " % locals()
+      
+      # If there are any NaN's in the results, issue a warning.
+      for mgNum, o in enumerate(results.xsilObjects):
+        for v in o.independentVariables:
+          if numpy.isnan(v['array']).any():
+            print >> sys.stderr, "Warning: Coordinate '%s' in moment group %i of file '%s' contains a NaN." % (v['name'], mgNum+1, sourceFile)
+        for v in o.dependentVariables:
+          if numpy.isnan(v['array']).any():
+            print >> sys.stderr, "Warning: Dependent variable '%s' in moment group %i of file '%s' contains a NaN." % (v['name'], mgNum+1, sourceFile)
+      
+      resultsFileContents = file(resultsFullPath).read()
+      
+      for xsilObject in results.xsilObjects:
+        if hasattr(xsilObject.data, 'filename'):
+          # If the moment group has a data file name, then we need to copy it to the expected results file
+          newDataFilename = xsilObject.data.filename.replace(os.path.splitext(sourceFile)[0], os.path.splitext(expectedResultsFile)[0], 1)
+          
+          resultsFileContents = resultsFileContents.replace(xsilObject.data.filename, newDataFilename)
+          
+          shutil.copyfile(os.path.join(testDir, xsilObject.data.filename),
+                          os.path.join(os.path.split(absPath)[0], newDataFilename))
+      
+      file(expectedResultsFullPath, 'w').write(resultsFileContents)
+    else:
+      expectedResults = XSILFile(expectedResultsFullPath)
+      
+      self.assert_(len(results.xsilObjects) == len(expectedResults.xsilObjects))
+      
+      momentGroupElements = xsilFileElement.getChildElementsByTagName('moment_group', optional=True)
+      if momentGroupElements:
+        self.assert_(len(momentGroupElements) == len(results.xsilObjects))
+      else:
+        momentGroupElements = [None]*len(results.xsilObjects)
+      
+      for mgNum, (o1, o2, mgElem) in enumerate(zip(results.xsilObjects, expectedResults.xsilObjects, momentGroupElements)):
+        currentAbsoluteTolerance = absoluteTolerance
+        currentRelativeTolerance = relativeTolerance
+        self.assert_(len(o1.independentVariables) == len(o2.independentVariables),
+                     "The number of independent variables in moment group %(mgNum)i doesn't match." % locals())
+        self.assert_(len(o1.dependentVariables) == len(o2.dependentVariables),
+                     "The number of dependent variables in moment group %(mgNum)i doesn't match." % locals())
+        
+        if mgElem:
+          if mgElem.hasAttribute('absolute_tolerance'):
+            currentAbsoluteTolerance = float(mgElem.getAttribute('absolute_tolerance'))
+          if mgElem.hasAttribute('relative_tolerance'):
+            currentRelativeTolerance = float(mgElem.getAttribute('relative_tolerance'))
+        
+        self.assert_(currentAbsoluteTolerance != None and currentRelativeTolerance != None, "An absolute and a relative tolerance must be specified.")
+        
+        for v1, v2 in zip(o1.independentVariables, o2.independentVariables):
+          self.assert_(v1['name'] == v2['name'])
+          self.assert_(v1['length'] == v2['length'])
+          # These are the coordinates, we just specify a constant absolute and relative tolerance.
+          # No-one should need to change these
+          self.assert_(array_approx_equal(v1['array'], v2['array'], 1e-7, 1e-6),
+                      "Coordinate '%s' in moment group %i of file '%s' didn't pass tolerance criteria." % (v1['name'], mgNum+1, sourceFile))
+        
+        for v1, v2 in zip(o1.dependentVariables, o2.dependentVariables):
+          self.assert_(v1['name'] == v2['name'])
+          self.assert_(pass_nan_test(v1['array'], v2['array']),
+                       "Dependent variable '%s' in moment group %i of file '%s' had a NaN where the expected results didn't (or vice-versa)." % (v1['name'], mgNum+1, sourceFile))
+          self.assert_(array_approx_equal(v1['array'], v2['array'], currentAbsoluteTolerance, currentRelativeTolerance),
+                       "Dependent variable '%s' in moment group %i of file '%s' failed to pass tolerance criteria." % (v1['name'], mgNum+1, sourceFile))
+  
+  # Test has succeeded, so save our checksum for the source file and copy the source file
+  file(checksumFilePath, 'w').write(currentChecksum)
+  
+  lastKnownGoodSourcePath = os.path.join(testDir, simulationName + '_last_known_good.cc')
+  file(lastKnownGoodSourcePath, 'w').write(sourceContents)
+
+def partial(func, *args, **keywords):
+  def newfunc(*fargs, **fkeywords):
+    newkeywords = keywords.copy()
+    newkeywords.update(fkeywords)
+    return func(*(args + fargs), **newkeywords)
+  return newfunc
+
+
+def main(argv=None):
+  if argv is None:
+    argv = sys.argv
+  try:
+    try:
+      opts, args = getopt.getopt(argv[1:], "ho:v", ["help", "output="])
+    except getopt.error, msg:
+      raise Usage(msg)
+  
+    # option processing
+    for option, value in opts:
+      if option == "-v":
+        verbose = True
+      if option in ("-h", "--help"):
+        raise Usage(help_message)
+      if option in ("-o", "--output"):
+        output = value
+  
+  except Usage, err:
+    print >> sys.stderr, sys.argv[0].split("/")[-1] + ": " + str(err.msg)
+    print >> sys.stderr, "\t for help use --help"
+    return 2
+  
+  basePath = os.path.dirname(__file__)
+  
+  resultsPath = os.path.join(basePath, 'testsuite_results')
+  if not os.path.exists(resultsPath):
+    os.mkdir(resultsPath)
+  resultsPath = os.path.abspath(resultsPath)
+  
+  print "Saving test results in %(resultsPath)s" % locals()
+  
+  testsuites = {}
+  baseSuiteName = 'testsuite'
+  baseSuitePath = os.path.join(basePath, baseSuiteName)
+  
+  for root, dirs, files in os.walk(baseSuitePath):
+    # First remove directories we don't want to traverse
+    for dirName in ['.svn']:
+      if dirName in dirs:
+        dirs.remove(dirName)
+    # Remove the 'testsuite/' part of the path
+    dirRelativeToBase = root[(len(baseSuitePath)+1):]
+    if dirRelativeToBase:
+      testSuiteName = os.path.join(baseSuiteName, dirRelativeToBase)
+    else:
+      testSuiteName = baseSuiteName
+    
+    # If we have .xmds files in this path, then create a TestCase subclass
+    xmdsTestScripts = [filename for filename in files if os.path.splitext(filename)[1].lower() == '.xmds']
+    
+    if xmdsTestScripts:
+      class ScriptTestCase(unittest.TestCase):
+        # Create test functions for each test script using 'scriptTestingFunction'
+        # These test function names are of the form 'test_ScriptName'
+        for scriptName in xmdsTestScripts:
+          prefix = os.path.splitext(scriptName)[0]
+          absPath = os.path.abspath(os.path.join(root, scriptName))
+          testDir = os.path.join(resultsPath, dirRelativeToBase)
+          locals()['test_' + prefix] = partial(scriptTestingFunction, root, scriptName, testDir, absPath)
+          locals()['test_' + prefix].__doc__ = os.path.join(dirRelativeToBase, scriptName)
+      
+      # Create a TestSuite from that class
+      suite = unittest.defaultTestLoader.loadTestsFromTestCase(ScriptTestCase)
+      testsuites[testSuiteName] = suite
+    
+    if not testSuiteName in testsuites:
+      testsuites[testSuiteName] = unittest.TestSuite()
+    
+    suite = testsuites[testSuiteName]
+    # Add our TestSuite as a sub-suite of all parent suites
+    head = testSuiteName
+    while True:
+      head, tail = os.path.split(head)
+      if not head or not tail:
+        break
+      testsuites[head].addTest(suite)
+  
+  
+  suitesToRun = list()
+  if len(args):
+    for suiteName in args:
+      fullSuiteName = os.path.join(baseSuiteName, suiteName)
+      if fullSuiteName in testsuites:
+        suitesToRun.append(testsuites[fullSuiteName])
+      else:
+        print >> sys.stderr, "Unable to find test '%(suiteName)s'" % locals()
+  else:
+    suitesToRun.append(testsuites[baseSuiteName])
+  suitesToRun.append(unittest.defaultTestLoader.loadTestsFromModule(CodeParser))
+  
+  fullSuite = unittest.TestSuite(tests=suitesToRun)
+  
+  return not unittest.TextTestRunner().run(fullSuite).wasSuccessful()
+
+
+if __name__ == "__main__":
+  sys.exit(main())
diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb.xmds b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb.xmds
new file mode 100644
index 0000000..8963ae5
--- /dev/null
+++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb.xmds
@@ -0,0 +1,667 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+	<testing>
+		<arguments> --Ndens=1e15 --Lcell=10.0e-2  --Temperature=1e-3  --Pwidth=0.4e-6  --delta1=0 --delta2=0 --delta3=0  --E1o=0 --E2o=1e2 --E3o=0 --E4o=0</arguments>
+		<xsil_file name="realistic_Rb.xsil" expected="realistic_Rb_expected.xsil" absolute_tolerance="1e-7" relative_tolerance="1e-5">
+			<moment_group number="0" absolute_tolerance="1e-7" relative_tolerance="1e-6" />
+		</xsil_file>
+	</testing>
+
+	<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="estimate" 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="2"   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[
+					// 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);
+
+					dE1_dz = 0.16666666666666666*i*eta1*(2.449489742783178*r0903 + 4.242640687119286*r1004 + 6.*r1105) - Lt[E1];
+					dE2_dz = -0.4082482904638631*i*eta1*(r0907 + r1008) - Lt[E2];
+					dE3_dz = -0.3333333333333333*i*eta2*(1.7320508075688772*r1406 + 3.*r1507 + 4.242640687119286*r1608) - Lt[E3];
+					dE4_dz = (i*eta2*(2.449489742783178*r1301 + 3*r1402 + 3*r1503 + 2.449489742783178*r1604))/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: 
+-->
diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil
new file mode 100644
index 0000000..409374d
--- /dev/null
+++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil
@@ -0,0 +1,699 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+        <testing>
+                <arguments> --Ndens=1e15 --Lcell=10.0e-2  --Temperature=1e-3  --Pwidth=0.4e-6  --delta1=0 --delta2=0 --delta3=0  --E1o=0 --E2o=1e2 --E3o=0 --E4o=0</arguments>
+                <xsil_file name="realistic_Rb.xsil" expected="realistic_Rb_expected.xsil" absolute_tolerance="1e-7" relative_tolerance="1e-5">
+                        <moment_group number="0" absolute_tolerance="1e-7" relative_tolerance="1e-6" />
+                </xsil_file>
+        </testing>
+
+        <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="estimate" 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="2"   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[
+                                        // 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);
+
+                                        dE1_dz = 0.16666666666666666*i*eta1*(2.449489742783178*r0903 + 4.242640687119286*r1004 + 6.*r1105) - Lt[E1];
+                                        dE2_dz = -0.4082482904638631*i*eta1*(r0907 + r1008) - Lt[E2];
+                                        dE3_dz = -0.3333333333333333*i*eta2*(1.7320508075688772*r1406 + 3.*r1507 + 4.242640687119286*r1608) - Lt[E3];
+                                        dE4_dz = (i*eta2*(2.449489742783178*r1301 + 3*r1402 + 3*r1503 + 2.449489742783178*r1604))/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>
+
+
+<info>
+Script compiled with XMDS2 version VERSION_PLACEHOLDER (SUBVERSION_REVISION_PLACEHOLDER)
+See http://www.xmds.org for more information.
+
+Variables that can be specified on the command line:
+  Command line argument E1o = 0.000000e+00
+  Command line argument E2o = 1.000000e+02
+  Command line argument E3o = 0.000000e+00
+  Command line argument E4o = 0.000000e+00
+  Command line argument delta1 = 0.000000e+00
+  Command line argument delta2 = 0.000000e+00
+  Command line argument delta3 = 0.000000e+00
+  Command line argument Pwidth = 4.000000e-07
+  Command line argument Lcell = 1.000000e-01
+  Command line argument Ndens = 1.000000e+15
+  Command line argument Temperature = 1.000000e-03
+</info>
+
+<XSIL Name="moment_group_1">
+  <Param Name="n_independent">2</Param>
+  <Array Name="variables" Type="Text">
+    <Dim>6</Dim>
+    <Stream><Metalink Format="Text" Delimiter=" \n"/>
+z t I1_out I2_out I3_out I4_out 
+    </Stream>
+  </Array>
+  <Array Name="data" Type="double">
+    <Dim>101</Dim>
+    <Dim>1000</Dim>
+    <Dim>6</Dim>
+    <Stream><Metalink Format="Binary" UnsignedLong="uint32" precision="double" Type="Remote" Encoding="LittleEndian"/>
+realistic_Rb_expected_mg0.dat
+    </Stream>
+  </Array>
+</XSIL>
+</simulation>
diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat
new file mode 100644
index 0000000..ff30ab5
--- /dev/null
+++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat