added code from newly regenerated mathematica->xmds files

* the main change is presense of the header and footer lines which
  allows to easily remove add this code at the right amount

1 files changed, 93 insertions, 113 deletions

diff --git a/xmds2/realistic_Rb/realistic_Rb.xmds b/xmds2/realistic_Rb/realistic_Rb.xmds
index 1111b33..94b398c 100644
--- a/xmds2/realistic_Rb/realistic_Rb.xmds
+++ b/xmds2/realistic_Rb/realistic_Rb.xmds
@@ -55,19 +55,21 @@
 
 				// 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;
+//---------------- Constants.txt starts ------------------
+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;
+//---------------- Constants.txt ends   ------------------
 
 				const double c=3.e8;
 				const double k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K]
@@ -225,38 +227,10 @@
 
 	<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
+			<!-- read from Mathematica generated RbChosenRho.txt -->
+<!-- ############### RbChosenRho.txt starts ################ -->
+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 
+<!-- ############### RbChosenRho.txt ends  ################# -->
 		</components>
 		<initialisation>
 			<!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
@@ -272,38 +246,40 @@
 			//       through population normalization
 
 			//read from Mathematica generated RbInits.txt
-				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;
+//---------------- RbInits.txt starts ------------------
+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;
+//---------------- RbInits.txt ends  ------------------
 			]]>
 		</initialisation>
 	</vector>
@@ -344,38 +320,40 @@
 					<![CDATA[
 						// Equations of motions according to Simon's mathematica code
 						//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));
+//---------------- RbEquations.txt starts ------------------
+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));
+//---------------- RbEquations.txt ends  ------------------
 					]]>
         </operator>
 				<!--
@@ -394,10 +372,12 @@
 				<dependencies>density_matrix</dependencies>
 				<![CDATA[
 				//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];
+//---------------- RbPropEquations.txt starts ------------------
+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];
+//---------------- RbPropEquations.txt ends  ------------------
 				]]>
 			</operators>
 		</integrate>