added new test case using --gt argument

2 files changed, 89 insertions, 76 deletions

diff --git a/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected.xsil b/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected.xsil
index c09460c..66b3827 100644
--- a/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected.xsil
+++ b/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected.xsil
@@ -1,9 +1,15 @@
-<?xml version="1.0"?>
-<simulation xmds-version="2">
+<?xml version="1.0" ?><!-- Below text defines necessary includes files generated by mathematica 
+                 and references to them --><!DOCTYPE simulation [
+<!ENTITY RbAtomConstantsFile SYSTEM "RbAtomConstants.cdata" >
+<!ENTITY RbChosenRhoFile SYSTEM "RbChosenRho.txt" >
+<!ENTITY RbInitsFile SYSTEM "RbInits.cdata" >
+<!ENTITY RbEquationsFile SYSTEM "RbEquations.cdata" >
+<!ENTITY RbPropEquationsFile SYSTEM "RbPropEquations.cdata" >
+]><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  --Ep1o=1e4 --Ep2o=1e2 --Ep3o=1e1 --Ep4o=1e1 --Em1o=1e4 --Em2o=1e2 --Em3o=1e1 --Em4o=1e1 --WLx=0 --WLy=0 --WLz=0</arguments>
-                <xsil_file name="realistic_Rb_and_fields.xsil" expected="realistic_Rb_and_fields_expected.xsil" absolute_tolerance="1e-7" relative_tolerance="1e-5">
-                        <moment_group number="0" absolute_tolerance="1e-7" relative_tolerance="1e-6" />
+                <arguments> --Ndens=1e15 --gt=0.1e6 --Lcell=10.0e-2  --Temperature=1e-3  --Pwidth=0.4e-6  --delta1=0 --delta2=0 --delta3=0  --Ep1o=1e4 --Ep2o=1e2 --Ep3o=1e1 --Ep4o=1e1 --Em1o=1e4 --Em2o=1e2 --Em3o=1e1 --Em4o=1e1 --WLx=0 --WLy=0 --WLz=0</arguments>
+                <xsil_file absolute_tolerance="1e-7" expected="realistic_Rb_and_fields_expected.xsil" name="realistic_Rb_and_fields.xsil" relative_tolerance="1e-5">
+                        <moment_group absolute_tolerance="1e-7" number="0" relative_tolerance="1e-6"/>
                 </xsil_file>
         </testing>
 
@@ -22,16 +28,16 @@
                 resonance as fields E2 and E1.
 
 
-                *              --------------- | F=1, 2P_3/2 >
+                *              --------------- | F=1, 2P_3/2 &gt;
                 *                \        \    
-                *                 \ E3_r   \    -------- | F=2, 2P_+1/2 >
+                *                 \ E3_r   \    -------- | F=2, 2P_+1/2 &gt;
                 *                  \   E4_r \   /    \
                 *                   \        \ / E2_l \ 
                 *                    \        /        \ E1_l
-                *   | F=2, 2S_1/2 > --------------      \
+                *   | F=2, 2S_1/2 &gt; --------------      \
                 *                               \        \
                 *                                \        \
-                *                               ------------- | F=1, 2S_1/2 >
+                *                               ------------- | F=1, 2S_1/2 &gt;
                 *            
 
 
@@ -49,13 +55,10 @@
 
         <features>
                 <globals>
+                        <!--atom related constants-->
+                        <!--read from Mathematica generated RbAtomConstants.cdata-->
                         <![CDATA[
-                                // Some numerical constants
-                                const double pi = M_PI; 
-
-                                // atom related constants
-                                //read from Mathematica generated Constants.txt
-//---------------- Constants.txt starts ------------------
+//---------------- RbAtomConstants.cdata starts ------------------
 const double ha0 = 2.1471788680034824e10;
 const double ha1 = 2.558764384495815e9;
 const double g1 = 3.612847284945266e7;
@@ -69,7 +72,13 @@ const double eta2 = 5.5397657647874e-6;
 const double rt6 = 2.449489742783178;
 const double rt3 = 1.7320508075688772;
 const double rt2 = 1.4142135623730951;
-//---------------- Constants.txt ends   ------------------
+//---------------- RbAtomConstants.cdata ends   ------------------
+]]>
+
+
+                        <![CDATA[
+                                // Some numerical constants
+                                const double pi = M_PI; 
 
                                 const double c=3.e8;
                                 const double k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K]
@@ -90,46 +99,44 @@ const double rt2 = 1.4142135623730951;
                                 // 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);
-
                                 // 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 />
+                <benchmark/>
     <arguments>
                         <!-- Rabi frequency divided by 2 in [1/s] -->
-      <argument name="Ep1o" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
-      <argument name="Em1o" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
-      <argument name="Ep2o" type="real" default_value="0.05*(2*M_PI*1e6)" />
-      <argument name="Em2o" type="real" default_value="0.05*(2*M_PI*1e6)" />
-      <argument name="Ep3o" type="real" default_value="2*3.0*(2*M_PI*1e6)" />
-      <argument name="Em3o" type="real" default_value="2*3.0*(2*M_PI*1e6)" />
-      <argument name="Ep4o" type="real" default_value=".01*(2*M_PI*1e6)" />
-      <argument name="Em4o" type="real" default_value=".01*(2*M_PI*1e6)" />
+      <argument default_value="2*1.5*(2*M_PI*1e6)" name="Ep1o" type="real"/>
+      <argument default_value="2*1.5*(2*M_PI*1e6)" name="Em1o" type="real"/>
+      <argument default_value="0.05*(2*M_PI*1e6)" name="Ep2o" type="real"/>
+      <argument default_value="0.05*(2*M_PI*1e6)" name="Em2o" type="real"/>
+      <argument default_value="2*3.0*(2*M_PI*1e6)" name="Ep3o" type="real"/>
+      <argument default_value="2*3.0*(2*M_PI*1e6)" name="Em3o" type="real"/>
+      <argument default_value=".01*(2*M_PI*1e6)" name="Ep4o" type="real"/>
+      <argument default_value=".01*(2*M_PI*1e6)" name="Em4o" type="real"/>
                         <!-- 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" />
+      <argument default_value="0.0" name="delta1" type="real"/>
+      <argument default_value="0.0" name="delta2" type="real"/>
+      <argument default_value="0.0" name="delta3" type="real"/>
                         <!--Pulse duration/width [s] -->
-                        <argument name="Pwidth"  type="real" default_value="0.1e-6" />
+                        <argument default_value="0.1e-6" name="Pwidth" type="real"/>
                         <!--  Atom and cell properties -->
+                        <!--repopulation rate (atoms flying in/out the laser beam)  in [1/s]-->
+                        <argument default_value="0.1e6" name="gt" type="real"/>
                         <!--Cell length [m] -->
-                        <argument name="Lcell"  type="real" default_value="1.5e-2" />
+                        <argument default_value="1.5e-2" name="Lcell" type="real"/>
                         <!--Density of atoms [1/m^3] -->
-                        <argument name="Ndens"  type="real" default_value="1e15" />
+                        <argument default_value="1e15" name="Ndens" type="real"/>
                         <!-- Projections of magnetic field in Larmor frequency units [1/s] -->
-                        <argument name="WLx"  type="real" default_value="0" />
-                        <argument name="WLy"  type="real" default_value="0" />
-                        <argument name="WLz"  type="real" default_value="0" />
+                        <argument default_value="0" name="WLx" type="real"/>
+                        <argument default_value="0" name="WLy" type="real"/>
+                        <argument default_value="0" name="WLz" type="real"/>
                         <!--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" />
+                        <argument default_value="5" name="Temperature" type="real"/>
                         <!-- 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[
@@ -142,14 +149,14 @@ const double rt2 = 1.4142135623730951;
                                 V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; 
                         ]]>
     </arguments>
-                <bing />
-                <diagnostics /> 
-                <fftw plan="estimate" threads="1" />
+                <bing/>
+                <diagnostics/> 
+                <fftw plan="estimate" threads="1"/>
                 <!--<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 />
+                <auto_vectorise/>
+                <halt_non_finite/>
         </features>
 
         <!-- 'z', 't', and 'v'  to have dimensions [m], [s], and [m/s]   -->
@@ -168,13 +175,13 @@ const double rt2 = 1.4142135623730951;
                                                 * 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)" />
+                        <dimension domain="(-1e-6, 1e-6)" lattice="10000" name="t"/>
+                        <dimension domain="(V_maxwell_min, V_maxwell_max)" lattice="2" name="v"/>
                 </transverse_dimensions>
         </geometry>
 
         <!-- Rabi frequency --> 
-        <vector name="E_field" type="complex" initial_space="t">
+        <vector initial_basis="t" initial_space="t" name="E_field" type="complex">
                  <!--plus/minus circular polarization components-->
                 <components>Ep1 Em1 Ep2 Em2 Ep3 Em3 Ep4 Em4</components>
                 <initialisation>
@@ -194,33 +201,33 @@ const double rt2 = 1.4142135623730951;
         </vector>
 
         <!--Maxwell distribution probability p(v)-->
-        <computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real">
+        <vector dimensions="v" name="Maxwell_distribution_probabilities" type="real">
                 <components>probability_v</components>
-                <evaluation>
+                <initialisation>
                         <![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>
+                </initialisation>
+        </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">
+        <vector dimensions="" name="Maxwell_distribution_probabilities_norm" type="real">
                 <components>probability_v_norm</components>
-                <evaluation>
+                <initialisation>
                         <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>
+                </initialisation>
+        </vector>
 
 
         <!-- Averaged across Maxwell distribution fields amplitudes -->
-        <computed_vector name="E_field_avgd" dimensions="t" type="complex">
+        <computed_vector dimensions="t" name="E_field_avgd" type="complex">
                 <components>Ep1a Em1a Ep2a Em2a Ep3a Em3a Ep4a Em4a</components>
                 <evaluation>
                         <dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
@@ -238,12 +245,13 @@ const double rt2 = 1.4142135623730951;
                 </evaluation>
         </computed_vector>
 
-        <vector name="density_matrix" type="complex" initial_space="t">
+        <vector initial_basis="t" initial_space="t" name="density_matrix" type="complex">
                 <components>
                         <!-- read from Mathematica generated RbChosenRho.txt -->
-<!-- ############### RbChosenRho.txt starts ################ -->
+                        <!-- ############### RbChosenRho.txt starts ################ -->
 r0101 r0109 r0110 r0111 r0112 r0113 r0114 r0115 r0116 r0201 r0202 r0209 r0210 r0211 r0212 r0213 r0214 r0215 r0216 r0301 r0302 r0303 r0309 r0310 r0311 r0312 r0313 r0314 r0315 r0316 r0401 r0402 r0403 r0404 r0409 r0410 r0411 r0412 r0413 r0414 r0415 r0416 r0501 r0502 r0503 r0504 r0505 r0509 r0510 r0511 r0512 r0513 r0514 r0515 r0516 r0601 r0602 r0603 r0604 r0605 r0606 r0609 r0610 r0611 r0612 r0613 r0614 r0615 r0616 r0701 r0702 r0703 r0704 r0705 r0706 r0707 r0709 r0710 r0711 r0712 r0713 r0714 r0715 r0716 r0801 r0802 r0803 r0804 r0805 r0806 r0807 r0808 r0809 r0810 r0811 r0812 r0813 r0814 r0815 r0816 r0909 r0912 r0913 r0914 r0915 r0916 r1009 r1010 r1012 r1013 r1014 r1015 r1016 r1109 r1110 r1111 r1112 r1113 r1114 r1115 r1116 r1212 r1312 r1313 r1412 r1413 r1414 r1512 r1513 r1514 r1515 r1612 r1613 r1614 r1615 r1616 
 <!-- ############### RbChosenRho.txt ends  ################# -->
+
                 </components>
                 <initialisation>
                         <!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
@@ -257,9 +265,10 @@ r0101 r0109 r0110 r0111 r0112 r0113 r0114 r0115 r0116 r0201 r0202 r0209 r0210 r0
                         // TODO: Fix above. Make the equation of motion for r11 
                         //       and express other level, let's say r44
                         //       through population normalization
-
-                        //read from Mathematica generated RbInits.txt
-//---------------- RbInits.txt starts ------------------
+                        ]]>
+                        <!--read from Mathematica generated RbInits.cdata-->
+                        <![CDATA[
+//---------------- RbInits.cdata starts ------------------
 r0101 = 0.125;
 r0109 = 0;
 r0110 = 0;
@@ -396,15 +405,16 @@ r1613 = 0;
 r1614 = 0;
 r1615 = 0;
 r1616 = 0;
-//---------------- RbInits.txt ends  ------------------
-                        ]]>
+//---------------- RbInits.cdata ends  ------------------
+]]>
+
                 </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="ARK45" interval="Lcell" tolerance="1e-5"> 
                 <!--<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 than RK4 with 1000 steps-->
@@ -419,7 +429,7 @@ r1616 = 0;
                         <!--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">
+        <operator algorithm="SI" kind="cross_propagation" propagation_dimension="t">
                                         <integration_vectors>density_matrix</integration_vectors>
           <dependencies>E_field_avgd</dependencies>
           <boundary_condition kind="left">
@@ -434,10 +444,10 @@ r1616 = 0;
             ]]>
                                                 -->
           </boundary_condition>
+                                        <!-- Equations of motions according to Simon's mathematica code -->
+                                        <!-- read from Mathematica generated RbEquations.cdata -->
                                         <![CDATA[
-                                                // Equations of motions according to Simon's mathematica code
-                                                //read from Mathematica generated RbEquations.txt
-//---------------- RbEquations.txt starts ------------------
+//---------------- RbEquations.cdata starts ------------------
 dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1212)/3. + (g2*r1313)/6. - i*((conj(Em1a)*conj(r0109))/4. - (conj(Em4a)*conj(r0113))/(4.*rt6) - (Em1a*r0109)/4. + (Em4a*r0113)/(4.*rt6) + r0201*(WLx/2. - (i*WLy)/2.) - conj(r0201)*(WLx/2. + (i*WLy)/2.));
 dr0109_dt = (-(gt*r0109) - (gt + g1)*r0109)/2. - i*(-(conj(Ep1a)*conj(r0301))/(4.*rt6) + (conj(Ep2a)*conj(r0701))/(4.*rt6) - (conj(Em4a)*conj(r0913))/(4.*rt6) - (conj(Em1a)*r0101)/4. + (conj(Em1a)*r0909)/4. + r0209*(WLx/2. - (i*WLy)/2.) - r0110*(-WLx/(6.*rt2) - (i*WLy)/(6.*rt2)) - r0109*(-delta1 - v*Kvec1 - WLz/6.) + r0109*WLz);
 dr0110_dt = (-(gt*r0110) - (gt + g1)*r0110)/2. - i*(-(conj(Em1a)*conj(r0201))/(4.*rt2) - (conj(Ep1a)*conj(r0401))/(4.*rt2) - (conj(Em2a)*conj(r0601))/(4.*rt6) + (conj(Ep2a)*conj(r0801))/(4.*rt6) + (conj(Em1a)*conj(r1009))/4. - (conj(Em4a)*conj(r1013))/(4.*rt6) + (delta1 + v*Kvec1)*r0110 + r0210*(WLx/2. - (i*WLy)/2.) - r0111*(-WLx/(6.*rt2) - (i*WLy)/(6.*rt2)) - r0109*(-WLx/(6.*rt2) + (i*WLy)/(6.*rt2)) + r0110*WLz);
@@ -574,8 +584,9 @@ dr1613_dt = -((gt + g2)*r1613) - i*((conj(Em4a)*conj(r0116))/(4.*rt6) - (conj(Ep
 dr1614_dt = -((gt + g2)*r1614) - i*((conj(Em4a)*conj(r0216))/8. - (conj(Ep4a)*conj(r0416))/8. - (conj(Em3a)*conj(r0616))/(8.*rt3) - (conj(Ep3a)*conj(r0816))/(8.*rt3) - (Em4a*r0414)/(4.*rt6) + (Em3a*r0814)/(4.*rt2) - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r1614 + r1514*((2*WLx)/3. + (2*i*WLy)/3.) - r1613*(sqrt(0.6666666666666666)*WLx - sqrt(0.6666666666666666)*i*WLy) - r1615*(sqrt(0.6666666666666666)*WLx + sqrt(0.6666666666666666)*i*WLy) + r1614*(-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3 - (4*WLz)/3.));
 dr1615_dt = -((gt + g2)*r1615) - i*((conj(Em4a)*conj(r0316))/8. - (conj(Ep4a)*conj(r0516))/(4.*rt6) - (conj(Em3a)*conj(r0716))/8. - (Em4a*r0415)/(4.*rt6) + (Em3a*r0815)/(4.*rt2) + r1515*((2*WLx)/3. + (2*i*WLy)/3.) - r1616*((2*WLx)/3. + (2*i*WLy)/3.) - r1614*(sqrt(0.6666666666666666)*WLx - sqrt(0.6666666666666666)*i*WLy) + r1615*(-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3 - (4*WLz)/3.) - r1615*(-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3 - (2*WLz)/3.));
 dr1616_dt = -((gt + g2)*r1616) - i*((conj(Em4a)*conj(r0416))/(4.*rt6) - (conj(Em3a)*conj(r0816))/(4.*rt2) - (Em4a*r0416)/(4.*rt6) + (Em3a*r0816)/(4.*rt2) - r1615*((2*WLx)/3. - (2*i*WLy)/3.) + conj(r1615)*((2*WLx)/3. + (2*i*WLy)/3.));
-//---------------- RbEquations.txt ends  ------------------
-                                        ]]>
+//---------------- RbEquations.cdata ends  ------------------
+]]>
+
         </operator>
                                 <!--
                                                          According to xmds2 docs operator kind="ip" should be faster
@@ -583,7 +594,7 @@ dr1616_dt = -((gt + g2)*r1616) - i*((conj(Em4a)*conj(r0416))/(4.*rt6) - (conj(Em
                                                          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 constant="yes" kind="ex" type="imaginary">
                                         <operator_names>Lt</operator_names>
                                         <![CDATA[
                                         Lt = -i/c*kt;
@@ -591,9 +602,9 @@ dr1616_dt = -((gt + g2)*r1616) - i*((conj(Em4a)*conj(r0416))/(4.*rt6) - (conj(Em
                                 </operator>
         <integration_vectors>E_field</integration_vectors>
                                 <dependencies>density_matrix</dependencies>
+                                <!-- read from Mathematica generated RbPropEquations.cdata-->
                                 <![CDATA[
-                                //read from Mathematica generated RbPropEquations.txt
-//---------------- RbPropEquations.txt starts ------------------
+//---------------- RbPropEquations.cdata starts ------------------
 dEp1_dz = -((i*Ndens*(rt2*conj(r0309) + rt6*(conj(r0410) + rt2*conj(r0511)))*eta1)/(rt2*rt6)) - Lt[Ep1];
 dEm1_dz = -((i*Ndens*(rt2*rt6*conj(r0109) + rt6*conj(r0210) + rt2*conj(r0311))*eta1)/(rt2*rt6)) - Lt[Em1];
 dEp2_dz = (i*Ndens*(conj(r0709) + conj(r0810))*eta1)/rt6 - Lt[Ep2];
@@ -602,8 +613,9 @@ dEp3_dz = -((i*Ndens*(rt2*rt3*conj(r0612) + rt3*conj(r0713) + conj(r0814))*eta2)
 dEm3_dz = -((i*Ndens*(conj(r0614) + rt3*(conj(r0715) + rt2*conj(r0816)))*eta2)/rt3) - Lt[Em3];
 dEp4_dz = -((i*Ndens*(rt2*conj(r0212) + rt3*conj(r0313) + rt3*conj(r0414) + rt2*conj(r0515))*eta2)/rt3) - Lt[Ep4];
 dEm4_dz = (i*Ndens*(rt2*conj(r0113) + rt3*conj(r0214) + rt3*conj(r0315) + rt2*conj(r0416))*eta2)/rt3 - Lt[Em4];
-//---------------- RbPropEquations.txt ends  ------------------
-                                ]]>
+//---------------- RbPropEquations.cdata ends  ------------------
+]]>
+
                         </operators>
                 </integrate>
         </sequence>
@@ -611,7 +623,7 @@ dEm4_dz = (i*Ndens*(rt2*conj(r0113) + rt3*conj(r0214) + rt3*conj(r0315) + rt2*co
 
 
         <!-- The output to generate -->
-        <output format="binary" filename="realistic_Rb_and_fields.xsil">
+        <output filename="realistic_Rb_and_fields.xsil" format="binary">
                 <group>
       <sampling basis="t(1000) " initial_sample="yes">
                                 <dependencies>E_field_avgd</dependencies>
@@ -648,6 +660,7 @@ Variables that can be specified on the command line:
   Command line argument delta2 = 0.000000e+00
   Command line argument delta3 = 0.000000e+00
   Command line argument Pwidth = 4.000000e-07
+  Command line argument gt = 1.000000e+05
   Command line argument Lcell = 1.000000e-01
   Command line argument Ndens = 1.000000e+15
   Command line argument WLx = 0.000000e+00
diff --git a/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected_mg0.dat b/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected_mg0.dat
index c301bc7..768f1ee 100644
--- a/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected_mg0.dat
+++ b/xmds2/realistic_Rb_and_fields/tests/testsuite/realistic_Rb_and_fields_expected_mg0.dat