First draft of realistic Rb which compiles and runs

1 files changed, 243 insertions, 112 deletions

diff --git a/xmds2/realistic_Rb/realistic_Rb.xmds b/xmds2/realistic_Rb/realistic_Rb.xmds
index 6a40191..598715f 100644
--- a/xmds2/realistic_Rb/realistic_Rb.xmds
+++ b/xmds2/realistic_Rb/realistic_Rb.xmds
@@ -44,14 +44,26 @@
 	<features>
 		<globals>
 			<![CDATA[
+				// Some numerical constants
 				const double pi = M_PI; 
+				// proportional to splitting ratios  sqrt(6) , sqrt(3), sqrt(2)
+				const double rt6 = 2.449489742783178;
+				const double rt3 = 1.7320508075688772;
+				const double rt2 = 1.4142135623730951;
+
+
 				const double c=3.e8;
 				const double k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K]
 				const double lambda=794.7e-9; //wavelength in m
-				const double Kvec = 2*M_PI/lambda;  // k-vector
+				// Fields k-vector
+				const double Kvec = 2*M_PI/lambda;  
+				// Simplified k-vectors
+				const double Kvec1 = Kvec, Kvec2=Kvec, Kvec3=Kvec;
+
 				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
 				// eta will be calculated in the <arguments> section
 				double eta = 0;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				double eta1=0, eta2=0, eta3=0;
 				
 				//  ---------  Atom and cell properties -------------------------
 				// range of Maxwell distribution atomic velocities
@@ -66,18 +78,44 @@
 
 				// 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 G3=3.0 *(2*M_PI*1e6);
-				const double G4=3.0 *(2*M_PI*1e6);
+				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;
 
-				// total decay of i-th level branching ratios. Rij branching of i-th level to j-th
-				const double R41=0.5, R42=0.5;
-				const double R31=0.5, R32=0.5;
 
 
 				complex  E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies
 
-				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
+				// 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
@@ -121,6 +159,10 @@
 
 				// 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 />
@@ -211,38 +253,119 @@
 
 	<!-- Averaged across Maxwell distribution density matrix components -->
 	<computed_vector name="density_matrix_averaged" dimensions="t" type="complex">
-		<components>r11a r22a r33a r12a r13a r14a r23a r24a r34a r44a</components>
+		<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;
-			r11a=r11*prob_v_normalized;
-			r22a=r22*prob_v_normalized;
-			r33a=r33*prob_v_normalized;
-			r12a=r12*prob_v_normalized;
-			r13a=r13*prob_v_normalized;
-			r14a=r14*prob_v_normalized;
-			r23a=r23*prob_v_normalized;
-			r24a=r24*prob_v_normalized;
-			r34a=r34*prob_v_normalized;
-			r44a=r44*prob_v_normalized;
+
+				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>r11 r22 r33 r44 r12 r13 r14 r23 r24 r34   r21 r31 r41 r32 r42 r43</components>-->
-		<!--<components>r11 r22 r33 r44 r12 r13 r14 r23 r24 r34</components>-->
-		<components>r11 r22 r33 r12 r13 r14 r23 r24 r34 r44</components>
-		<!--
-				 note one of the level population is redundant since
-				 r11+r22+r33+r44=1
-				 so r11 is missing
-		-->
+		<components>
+				r0101
+				r0113
+				r0202
+				r0214
+				r0303
+				r0309
+				r0315
+				r0404
+				r0410
+				r0416
+				r0505
+				r0511
+				r0602
+				r0606
+				r0614
+				r0703
+				r0707
+				r0709
+				r0715
+				r0804
+				r0808
+				r0810
+				r0816
+				r0909
+				r0915
+				r1010
+				r1016
+				r1111
+				r1313
+				r1414
+				r1515
+				r1616
+		</components>
 		<initialisation>
 			<!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
-			<!-- Comment out no light field initial conditions 
 			<![CDATA[
 			// Note: 
 			// convergence is really slow if all populations concentrated at the bottom level |1>
@@ -253,67 +376,39 @@
 			// TODO: Fix above. Make the equation of motion for r11 
 			//       and express other level, let's say r44
 			//       through population normalization
-			r11 = 1; 
-			          r22 = 0; r33 = 0; r44 = 0;
-			 r12 = 0; r13 = 0; r14 = 0;
-			 r23 = 0; r24 = 0;
-			 r34 = 0;
+				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;
 			]]>
-			-->
-			<!-- Below initialization assumes strong E1 and E3 which were shining for long time before
-					 we even start to look at the problem -->
-			<!-- Precalculated by Nate via Mathematica steady state solver -->
-				<dependencies>E_field_avgd</dependencies>
-				<![CDATA[
-						E1ac = conj(E1a);
-						E2ac = conj(E2a);
-						E3ac = conj(E3a);
-						E4ac = conj(E4a);
-
-				// IMPORTANT: assumes no detunings
-				r11 = (gt*(4*mod2(E1a) + (G3 + gt)*(G3 + 2*gt))*((G4 +
-				2*gt)*(4*mod2(E3a) + gt*(G4 + gt)) + 4*mod2(E3a)*G4*(R41 -
-				R42)))/(2.*(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
-				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
-
-				r13 = -((E1ac*gt*(G3 + gt)*i*((G4 + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt))
-				+ 4*mod2(E3a)*G4*(R41 -
-				R42)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
-				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
-
-				r22 = (gt*(4*mod2(E3a) + (G4 + gt)*(G4 + 2*gt))*((G3 +
-				2*gt)*(4*mod2(E1a) + gt*(G3 + gt)) + 4*mod2(E1a)*G3*(-R31 +
-				R32)))/(2.*(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
-				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
-
-				r24 = -((E3ac*gt*(G4 + gt)*i*((G3 + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt))
-				+ 4*mod2(E1a)*G3*(-R31 +
-				R32)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
-				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42)));
-
-				r33 = (2*mod2(E1a)*gt*((G4 + 2*gt)*(4*mod2(E3a) + gt*(G4 + gt)) +
-				4*mod2(E3a)*G4*(R41 -
-				R42)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
-				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42));
-
-				r44 = (2*mod2(E3a)*gt*((G3 + 2*gt)*(4*mod2(E1a) + gt*(G3 + gt)) +
-				4*mod2(E1a)*G3*(-R31 +
-				R32)))/(-16*mod2(E1a)*mod2(E3a)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1a) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3a) +
-				gt*(G4 + gt)) - 4*mod2(E3a)*G4*R42) + 4*mod2(E1a)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3a) + gt*(G4 + gt))) + 4*mod2(E3a)*G4*R42));
-
-				]]>
 		</initialisation>
 	</vector>
 
@@ -330,7 +425,8 @@
 
 		<!--<integrate algorithm="SIC" interval="4e-2" steps="200">-->
 
-			<samples>100 100</samples>
+			<samples>100</samples>
+			<!--<samples>100 100</samples>-->
 			<!--Use the next line for debuging to see velocity dependence. Uncomment/switch on output groups 3,4-->
 			<!--<samples>100 100 100 100</samples>--> 
 			<operators>
@@ -355,26 +451,59 @@
 						E3ac = conj(E3a);
 						E4ac = conj(E4a);
 
-						r21=conj(r12);
-						r31=conj(r13);
-						r41=conj(r14);
-						r32=conj(r23);
-						r42=conj(r24);
-						r43=conj(r34);
-						//r44=1- r11 - r22 - r33;
-
-						// Equations of motions according to Nate's mathematica code
-						dr11_dt = gt/2 - gt*r11 + i*(-(E1a*r13) - E4a*r14 + E1ac*r31 + E4ac*r41) + G3*r33*R31 + G4*r44*R41;
-						dr12_dt = -(gt*r12) + i*((-(delta1 + Kvec*v) + (delta2 + Kvec*v))*r12 - E2a*r13 - E3a*r14 + E1ac*r32 + E4ac*r42);
-						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1ac*r11) - E2ac*r12 - (delta1 + Kvec*v)*r13 + E1ac*r33 + E4ac*r43);
-						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4ac*r11) - E3ac*r12 - ((delta1 + Kvec*v) - (delta2 + Kvec*v) + (delta3 + Kvec*v))*r14 + E1ac*r34 + E4ac*r44);
-						dr22_dt = gt/2 - gt*r22 + i*(-(E2a*r23) - E3a*r24 + E2ac*r32 + E3ac*r42) + G3*r33*R32 + G4*r44*R42;
-						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1ac*r21) - E2ac*r22 - (delta2 + Kvec*v)*r23 + E2ac*r33 + E3ac*r43);
-						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4ac*r21) - E3ac*r22 - (delta3 + Kvec*v)*r24 + E2ac*r34 + E3ac*r44);
-						dr33_dt = i*(E1a*r13 + E2a*r23 - E1ac*r31 - E2ac*r32) - (G3 + gt)*r33;
-
-						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1a*r14 + E2a*r24 - E4ac*r31 - E3ac*r32 + ((delta2 + Kvec*v) - (delta3 + Kvec*v))*r34);
-						dr44_dt = i*(E4a*r14 + E3a*r24 - E4ac*r41 - E3ac*r42) - (G4 + gt)*r44;
+						// Density matrix is Hermitian so we use r_ij=conj(r_ji)
+
+						r1301 = conj(r0113);
+						r1402 = conj(r0214);
+						r0903 = conj(r0309);
+						r1503 = conj(r0315);
+						r1004 = conj(r0410);
+						r1604 = conj(r0416);
+						r1105 = conj(r0511);
+						r0206 = conj(r0602);
+						r1406 = conj(r0614);
+						r0307 = conj(r0703);
+						r0907 = conj(r0709);
+						r1507 = conj(r0715);
+						r0408 = conj(r0804);
+						r1008 = conj(r0810);
+						r1608 = conj(r0816);
+						r1509 = conj(r0915);
+						r1610 = conj(r1016);
+
+						// Equations of motions according to Simon's mathematica code
+						dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1313)/6. - i*((r0113*E4a)/(4.*rt6) - (r1301*E4ac)/(4.*rt6));
+						dr0113_dt = (-(gt*r0113) - (gt + g2)*r0113)/2. - i*(WL*r0113 - ((2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0113 + (r0101*E4ac)/(4.*rt6) - (r1313*E4ac)/(4.*rt6));
+						dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (r1402*E4ac)/8.);
+						dr0214_dt = (-(gt*r0214) - (gt + g2)*r0214)/2. - i*((WL*r0214)/2. - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0214 - (r0206*E3ac)/(8.*rt3) + (r0202*E4ac)/8. - (r1414*E4ac)/8.);
+						dr0303_dt = gt/8. - gt*r0303 + (g1*r0909)/12. + (g1*r1010)/3. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1515)/4. - i*((r0309*E1a)/(4.*rt6) + (r0315*E4a)/8. - (r0903*E1ac)/(4.*rt6) - (r1503*E4ac)/8.);
+						dr0309_dt = (-(gt*r0309) - (gt + g1)*r0309)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0309) + (r0303*E1ac)/(4.*rt6) - (r0909*E1ac)/(4.*rt6) - (r0307*E2ac)/(4.*rt6) - (r1509*E4ac)/8.);
+						dr0315_dt = (-(gt*r0315) - (gt + g2)*r0315)/2. - i*(-(((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0315) - (r0915*E1ac)/(4.*rt6) - (r0307*E3ac)/8. + (r0303*E4ac)/8. - (r1515*E4ac)/8.);
+						dr0404_dt = gt/8. - gt*r0404 + (g1*r1010)/4. + (g1*r1111)/4. + (g2*r1414)/4. + (g2*r1515)/12. + (g2*r1616)/6. - i*((r0410*E1a)/(4.*rt2) + (r0416*E4a)/(4.*rt6) - (r1004*E1ac)/(4.*rt2) - (r1604*E4ac)/(4.*rt6));
+						dr0410_dt = (-(gt*r0410) - (gt + g1)*r0410)/2. - i*(-(WL*r0410)/2. + (delta1 + v*Kvec1)*r0410 + (r0404*E1ac)/(4.*rt2) - (r1010*E1ac)/(4.*rt2) - (r0408*E2ac)/(4.*rt6) - (r1610*E4ac)/(4.*rt6));
+						dr0416_dt = (-(gt*r0416) - (gt + g2)*r0416)/2. - i*(-(WL*r0416)/2. - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0416 - (r1016*E1ac)/(4.*rt2) - (r0408*E3ac)/(4.*rt2) + (r0404*E4ac)/(4.*rt6) - (r1616*E4ac)/(4.*rt6));
+						dr0505_dt = gt/8. - gt*r0505 + (g1*r1111)/2. + (g2*r1515)/6. + (g2*r1616)/3. - i*((r0511*E1a)/4. - (r1105*E1ac)/4.);
+						dr0511_dt = (-(gt*r0511) - (gt + g1)*r0511)/2. - i*(-(WL*r0511) - (WL/6. - delta1 - v*Kvec1)*r0511 + (r0505*E1ac)/4. - (r1111*E1ac)/4.);
+						dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (r1402*E3ac)/(8.*rt3));
+						dr0606_dt = gt/8. - gt*r0606 + (g1*r0909)/12. + (g1*r1010)/12. + (g2*r1313)/4. + (g2*r1414)/12. - i*(-(r0614*E3a)/(8.*rt3) + (r1406*E3ac)/(8.*rt3));
+						dr0614_dt = (-(gt*r0614) - (gt + g2)*r0614)/2. - i*((-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0614 - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0614 - (r0606*E3ac)/(8.*rt3) + (r1414*E3ac)/(8.*rt3) + (r0602*E4ac)/8.);
+						dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (r0903*E2ac)/(4.*rt6) + (r1503*E3ac)/8.);
+						dr0707_dt = gt/8. - gt*r0707 + (g1*r0909)/12. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1414)/3. + (g2*r1515)/4. - i*(-(r0709*E2a)/(4.*rt6) - (r0715*E3a)/8. + (r0907*E2ac)/(4.*rt6) + (r1507*E3ac)/8.);
+						dr0709_dt = (-(gt*r0709) - (gt + g1)*r0709)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0709) + (-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0709 + (r0703*E1ac)/(4.*rt6) - (r0707*E2ac)/(4.*rt6) + (r0909*E2ac)/(4.*rt6) + (r1509*E3ac)/8.);
+						dr0715_dt = (-(gt*r0715) - (gt + g2)*r0715)/2. - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0715 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0715 + (r0915*E2ac)/(4.*rt6) - (r0707*E3ac)/8. + (r1515*E3ac)/8. + (r0703*E4ac)/8.);
+						dr0804_dt = -(gt*r0804) - i*((WL*r0804)/2. + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0804 + (r0810*E1a)/(4.*rt2) + (r0816*E4a)/(4.*rt6) + (r1004*E2ac)/(4.*rt6) + (r1604*E3ac)/(4.*rt2));
+						dr0808_dt = gt/8. - gt*r0808 + (g1*r1010)/12. + (g1*r1111)/12. + (g2*r1414)/12. + (g2*r1515)/4. + (g2*r1616)/2. - i*(-(r0810*E2a)/(4.*rt6) - (r0816*E3a)/(4.*rt2) + (r1008*E2ac)/(4.*rt6) + (r1608*E3ac)/(4.*rt2));
+						dr0810_dt = (-(gt*r0810) - (gt + g1)*r0810)/2. - i*((delta1 + v*Kvec1)*r0810 + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0810 + (r0804*E1ac)/(4.*rt2) - (r0808*E2ac)/(4.*rt6) + (r1010*E2ac)/(4.*rt6) + (r1610*E3ac)/(4.*rt2));
+						dr0816_dt = (-(gt*r0816) - (gt + g2)*r0816)/2. - i*((WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0816 - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0816 + (r1016*E2ac)/(4.*rt6) - (r0808*E3ac)/(4.*rt2) + (r1616*E3ac)/(4.*rt2) + (r0804*E4ac)/(4.*rt6));
+						dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (r0903*E1ac)/(4.*rt6) - (r0907*E2ac)/(4.*rt6));
+						dr0915_dt = (-((gt + g1)*r0915) - (gt + g2)*r0915)/2. - i*((-WL/6. - delta1 - v*Kvec1)*r0915 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0915 - (r0315*E1a)/(4.*rt6) + (r0715*E2a)/(4.*rt6) - (r0907*E3ac)/8. + (r0903*E4ac)/8.);
+						dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (r1004*E1ac)/(4.*rt2) - (r1008*E2ac)/(4.*rt6));
+						dr1016_dt = (-((gt + g1)*r1016) - (gt + g2)*r1016)/2. - i*(-((delta1 + v*Kvec1)*r1016) - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r1016 - (r0416*E1a)/(4.*rt2) + (r0816*E2a)/(4.*rt6) - (r1008*E3ac)/(4.*rt2) + (r1004*E4ac)/(4.*rt6));
+						dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (r1105*E1ac)/4.);
+						dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (r1301*E4ac)/(4.*rt6));
+						dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (r1406*E3ac)/(8.*rt3) + (r1402*E4ac)/8.);
+						dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (r1507*E3ac)/8. + (r1503*E4ac)/8.);
+						dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (r1608*E3ac)/(4.*rt2) + (r1604*E4ac)/(4.*rt6));
 					]]>
         </operator>
 				<!--
@@ -392,10 +521,10 @@
         <integration_vectors>E_field</integration_vectors>
 				<dependencies>density_matrix</dependencies>
           <![CDATA[
-					dE1_dz = i*eta*conj(r13) +Lt[E1] ;
-					dE2_dz = i*eta*conj(r23) +Lt[E2] ;
-					dE3_dz = i*eta*conj(r24) +Lt[E3] ;
-					dE4_dz = i*eta*conj(r14) +Lt[E4] ;
+					dE1_dz = 0.16666666666666666*eta1*(2.449489742783178*r0309 + 4.242640687119286*r0410 + 6.*r0511) - Lt[E1];
+					dE2_dz = -0.8164965809277261*eta1*(r0709 + r0810) - Lt[E2];
+					dE3_dz = -1.*eta2*(1.7320508075688772*r0614 + 3.*r0715 + 4.242640687119286*r0816) - Lt[E3];
+					dE4_dz = (4*eta2*(2.449489742783178*r0113 + 3*r0214 + 3*r0315 + 2.449489742783178*r0416))/3. - Lt[E4];
           ]]>
 			</operators>
 		</integrate>
@@ -418,6 +547,7 @@
 			</sampling>
 		</group>
 
+		<!--
 		<group>
       <sampling basis="t(100) v(10)" initial_sample="yes">
 				<dependencies>density_matrix_averaged</dependencies>
@@ -449,6 +579,7 @@
 				]]>
 			</sampling>
 		</group>
+		-->
 
 		<!-- use the following two groups only for debuging 
 				 otherwise they are quite useless and have to much information