From 249431a9209a2b744d43cad8fbdf48f891c8a1f9 Mon Sep 17 00:00:00 2001
From: Eugeniy Mikhailov <evgmik@gmail.com>
Date: Tue, 25 Sep 2012 11:00:33 -0400
Subject: heavy use of Hermitian properties in dRho/dt calculations

* also read constants from mathematica generated files
---
 xmds2/realistic_Rb/realistic_Rb.xmds | 136 +++++++++++++----------------------
 1 file changed, 50 insertions(+), 86 deletions(-)

diff --git a/xmds2/realistic_Rb/realistic_Rb.xmds b/xmds2/realistic_Rb/realistic_Rb.xmds
index f50d26c..79f5a04 100644
--- a/xmds2/realistic_Rb/realistic_Rb.xmds
+++ b/xmds2/realistic_Rb/realistic_Rb.xmds
@@ -46,12 +46,21 @@
 			<![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 = 0.00007949788511562656;
+				const double lambda2 = 0.00007802412096860508;
 				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
@@ -79,16 +88,6 @@
 				// 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;
 
@@ -96,26 +95,6 @@
 
 				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
@@ -167,9 +146,10 @@
     </arguments>
 		<bing />
 		<diagnostics /> 
-		<fftw plan="patient" threads="1" />
+		<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 />
+		<!--<openmp />-->
 		<auto_vectorise />
 		<halt_non_finite />
 	</features>
@@ -191,7 +171,7 @@
 							works for Doppler averaging in up to 400K for Rb when lasers are zero detuned
 			 -->
 			<dimension name="t"   lattice="10000"   domain="(-1e-6, 1e-6)" />
-			<dimension name="v"   lattice="100"   domain="(V_maxwell_min, V_maxwell_max)" />
+			<dimension name="v"   lattice="2"   domain="(V_maxwell_min, V_maxwell_max)" />
 		</transverse_dimensions>
 	</geometry>
 
@@ -376,6 +356,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;
@@ -451,59 +433,40 @@
 						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));
+						//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)*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.);
+						dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (conj(r0214)*E4ac)/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)*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. - (conj(r0309)*E1ac)/(4.*rt6) - (conj(r0315)*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) - (conj(r0703)*E2ac)/(4.*rt6) - (conj(r0915)*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) - (conj(r0703)*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) - (conj(r0410)*E1ac)/(4.*rt2) - (conj(r0416)*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) - (conj(r0804)*E2ac)/(4.*rt6) - (conj(r1016)*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) - (conj(r0804)*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. - (conj(r0511)*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));
+						dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (conj(r0214)*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) + (conj(r0614)*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.);
+						dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (conj(r0309)*E2ac)/(4.*rt6) + (conj(r0315)*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. + (conj(r0709)*E2ac)/(4.*rt6) + (conj(r0715)*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) + (conj(r0915)*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));
+						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)*E2ac)/(4.*rt6) + (conj(r0416)*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) + (conj(r0810)*E2ac)/(4.*rt6) + (conj(r0816)*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) + (conj(r1016)*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));
+						dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (conj(r0309)*E1ac)/(4.*rt6) - (conj(r0709)*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) - (conj(r0709)*E3ac)/8. + (conj(r0309)*E4ac)/8.);
+						dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (conj(r0410)*E1ac)/(4.*rt2) - (conj(r0810)*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) - (conj(r0810)*E3ac)/(4.*rt2) + (conj(r0410)*E4ac)/(4.*rt6));
+						dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (conj(r0511)*E1ac)/4.);
+						dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (conj(r0113)*E4ac)/(4.*rt6));
+						dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (conj(r0614)*E3ac)/(8.*rt3) + (conj(r0214)*E4ac)/8.);
+						dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (conj(r0715)*E3ac)/8. + (conj(r0315)*E4ac)/8.);
+						dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (conj(r0816)*E3ac)/(4.*rt2) + (conj(r0416)*E4ac)/(4.*rt6));
 					]]>
         </operator>
 				<!--
@@ -521,10 +484,11 @@
         <integration_vectors>E_field</integration_vectors>
 				<dependencies>density_matrix</dependencies>
 				<![CDATA[
-					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*(2.449489742783178*conj(r0309) + 4.242640687119286*conj(r0410) + 6.*conj(r0511))*eta1 - Lt[E1];
+				dE2_dz = -0.4082482904638631*i*(conj(r0709) + conj(r0810))*eta1 - Lt[E2];
+				dE3_dz = -0.3333333333333333*i*(1.7320508075688772*conj(r0614) + 3.*conj(r0715) + 4.242640687119286*conj(r0816))*eta2 - Lt[E3];
+				dE4_dz = (i*(2.449489742783178*conj(r0113) + 3*conj(r0214) + 3*conj(r0315) + 2.449489742783178*conj(r0416))*eta2)/3. - Lt[E4];
 				]]>
 			</operators>
 		</integrate>
-- 
cgit v1.2.3