updated test cases since we use new etas separated for each transition

2 files changed, 56 insertions, 242 deletions

diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil
index 409374d..36f5939 100644
--- a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil
+++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil
@@ -36,7 +36,7 @@
 
 
                 We are solving 
-                        dE/dz+(1/c)*dE/dt=i*eta*rho_ij,  where j level is higher then i.
+                        dE/dz+(1/c)*dE/dt=i*eta*rho_ij,  where j level is higher than 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
@@ -52,12 +52,23 @@
                         <![CDATA[
                                 // Some numerical constants
                                 const double pi = M_PI; 
-                                // proportional to splitting ratios  sqrt(6) , sqrt(3), sqrt(2)
+
+                                // atom related constants
+                                //read from Mathematica generated Constants.txt
+                                const double ha0 = 2.1471788680034824e10;
+                                const double ha1 = 2.558764384495815e9;
+                                const double g1 = 3.612847284945266e7;
+                                const double ha2 = 5.323020344462938e8;
+                                const double hb2 = 7.85178251911697e7;
+                                const double g2 = 3.8117309832741246e7;
+                                const double lambda1 = 7.949788511562656e-7;
+                                const double lambda2 = 7.802412096860509e-7;
+                                const double eta1 = 5.450949336831401e-6;
+                                const double eta2 = 5.5397657647874e-6;
                                 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
@@ -66,11 +77,6 @@
                                 // 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] 
@@ -85,44 +91,10 @@
                                 // 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
 
@@ -162,18 +134,12 @@
                                 // 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" />
+                <!--<fftw plan="patient" threads="1" />-->
                 <!-- I don't see any speed up on 6 core CPU even if use threads="6" -->
                 <!--<openmp />-->
                 <auto_vectorise />
@@ -257,84 +223,6 @@
                 </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
@@ -382,6 +270,8 @@
                         // 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
                                 r0101 = 0.125;
                                 r0113 = 0;
                                 r0202 = 0.125;
@@ -424,7 +314,7 @@
                 <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-->
+                <!--When ARK45 works it about 3 times faster than 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-->
@@ -452,64 +342,40 @@
                                                 -->
           </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));
+                                                //read from Mathematica generated RbEquations.txt
+                                                dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1313)/6. - i*((r0113*E4a)/(4.*rt6) - (conj(r0113)*conj(E4a))/(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*conj(E4a))/(4.*rt6) - (r1313*conj(E4a))/(4.*rt6));
+                                                dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (conj(r0214)*conj(E4a))/8.);
+                                                dr0214_dt = (-(gt*r0214) - (gt + g2)*r0214)/2. - i*((WL*r0214)/2. - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0214 - (conj(r0602)*conj(E3a))/(8.*rt3) + (r0202*conj(E4a))/8. - (r1414*conj(E4a))/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. - (conj(r0309)*conj(E1a))/(4.*rt6) - (conj(r0315)*conj(E4a))/8.);
+                                                dr0309_dt = (-(gt*r0309) - (gt + g1)*r0309)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0309) + (r0303*conj(E1a))/(4.*rt6) - (r0909*conj(E1a))/(4.*rt6) - (conj(r0703)*conj(E2a))/(4.*rt6) - (conj(r0915)*conj(E4a))/8.);
+                                                dr0315_dt = (-(gt*r0315) - (gt + g2)*r0315)/2. - i*(-(((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0315) - (r0915*conj(E1a))/(4.*rt6) - (conj(r0703)*conj(E3a))/8. + (r0303*conj(E4a))/8. - (r1515*conj(E4a))/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) - (conj(r0410)*conj(E1a))/(4.*rt2) - (conj(r0416)*conj(E4a))/(4.*rt6));
+                                                dr0410_dt = (-(gt*r0410) - (gt + g1)*r0410)/2. - i*(-(WL*r0410)/2. + (delta1 + v*Kvec1)*r0410 + (r0404*conj(E1a))/(4.*rt2) - (r1010*conj(E1a))/(4.*rt2) - (conj(r0804)*conj(E2a))/(4.*rt6) - (conj(r1016)*conj(E4a))/(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*conj(E1a))/(4.*rt2) - (conj(r0804)*conj(E3a))/(4.*rt2) + (r0404*conj(E4a))/(4.*rt6) - (r1616*conj(E4a))/(4.*rt6));
+                                                dr0505_dt = gt/8. - gt*r0505 + (g1*r1111)/2. + (g2*r1515)/6. + (g2*r1616)/3. - i*((r0511*E1a)/4. - (conj(r0511)*conj(E1a))/4.);
+                                                dr0511_dt = (-(gt*r0511) - (gt + g1)*r0511)/2. - i*(-(WL*r0511) - (WL/6. - delta1 - v*Kvec1)*r0511 + (r0505*conj(E1a))/4. - (r1111*conj(E1a))/4.);
+                                                dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (conj(r0214)*conj(E3a))/(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) + (conj(r0614)*conj(E3a))/(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*conj(E3a))/(8.*rt3) + (r1414*conj(E3a))/(8.*rt3) + (r0602*conj(E4a))/8.);
+                                                dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (conj(r0309)*conj(E2a))/(4.*rt6) + (conj(r0315)*conj(E3a))/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. + (conj(r0709)*conj(E2a))/(4.*rt6) + (conj(r0715)*conj(E3a))/8.);
+                                                dr0709_dt = (-(gt*r0709) - (gt + g1)*r0709)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0709) + (-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0709 + (r0703*conj(E1a))/(4.*rt6) - (r0707*conj(E2a))/(4.*rt6) + (r0909*conj(E2a))/(4.*rt6) + (conj(r0915)*conj(E3a))/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*conj(E2a))/(4.*rt6) - (r0707*conj(E3a))/8. + (r1515*conj(E3a))/8. + (r0703*conj(E4a))/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) + (conj(r0410)*conj(E2a))/(4.*rt6) + (conj(r0416)*conj(E3a))/(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) + (conj(r0810)*conj(E2a))/(4.*rt6) + (conj(r0816)*conj(E3a))/(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*conj(E1a))/(4.*rt2) - (r0808*conj(E2a))/(4.*rt6) + (r1010*conj(E2a))/(4.*rt6) + (conj(r1016)*conj(E3a))/(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*conj(E2a))/(4.*rt6) - (r0808*conj(E3a))/(4.*rt2) + (r1616*conj(E3a))/(4.*rt2) + (r0804*conj(E4a))/(4.*rt6));
+                                                dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (conj(r0309)*conj(E1a))/(4.*rt6) - (conj(r0709)*conj(E2a))/(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) - (conj(r0709)*conj(E3a))/8. + (conj(r0309)*conj(E4a))/8.);
+                                                dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (conj(r0410)*conj(E1a))/(4.*rt2) - (conj(r0810)*conj(E2a))/(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) - (conj(r0810)*conj(E3a))/(4.*rt2) + (conj(r0410)*conj(E4a))/(4.*rt6));
+                                                dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (conj(r0511)*conj(E1a))/4.);
+                                                dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (conj(r0113)*conj(E4a))/(4.*rt6));
+                                                dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (conj(r0614)*conj(E3a))/(8.*rt3) + (conj(r0214)*conj(E4a))/8.);
+                                                dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (conj(r0715)*conj(E3a))/8. + (conj(r0315)*conj(E4a))/8.);
+                                                dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (conj(r0816)*conj(E3a))/(4.*rt2) + (conj(r0416)*conj(E4a))/(4.*rt6));
                                         ]]>
         </operator>
                                 <!--
@@ -527,29 +393,11 @@
         <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];
+                                //read from Mathematica generated RbPropEquations.txt
+                                dE1_dz = 0.16666666666666666*i*Ndens*(2.449489742783178*conj(r0309) + 4.242640687119286*conj(r0410) + 6.*conj(r0511))*eta1 - Lt[E1];
+                                dE2_dz = -0.4082482904638631*i*Ndens*(conj(r0709) + conj(r0810))*eta1 - Lt[E2];
+                                dE3_dz = -0.3333333333333333*i*Ndens*(1.7320508075688772*conj(r0614) + 3.*conj(r0715) + 4.242640687119286*conj(r0816))*eta2 - Lt[E3];
+                                dE4_dz = (i*Ndens*(2.449489742783178*conj(r0113) + 3*conj(r0214) + 3*conj(r0315) + 2.449489742783178*conj(r0416))*eta2)/3. - Lt[E4];
                                 ]]>
                         </operators>
                 </integrate>
@@ -572,40 +420,6 @@
                         </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) -->
diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat
index ff30ab5..f17a063 100644
--- a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat
+++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat