added complex electromagnetic field propagation, rho44 is differential equations now

1 files changed, 14 insertions, 20 deletions

diff --git a/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds b/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds
index 78ab368..f5d7dfa 100644
--- a/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds
+++ b/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds
@@ -117,7 +117,7 @@
 	<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</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
@@ -134,7 +134,7 @@
 			// 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;
+			r11 = 1; r22 = 0; r33 = 0; r44 = 0;
 			r12 = 0; r13 = 0; r14 = 0;
 			r23 = 0; r24 = 0;
 			r34 = 0;
@@ -156,7 +156,7 @@
           <dependencies>E_field</dependencies>
           <boundary_condition kind="left">
             <![CDATA[
-							r11 = 1; r22 = 0; r33 = 0;
+							r11 = 1; r22 = 0; r33 = 0; r44 = 0;
 							r12 = 0; r13 = 0; r14 = 0;
 							r23 = 0; r24 = 0;
 							r34 = 0;
@@ -173,25 +173,19 @@
 						r32=conj(r23);
 						r42=conj(r24);
 						r43=conj(r34);
-						r44=1- r11 - r22 - r33;
+						//r44=1- r11 - r22 - r33;
 
 						// Equations of motions according to Simon's mathematica code
-						dr11_dt = gt/2. - gt*r11 + E1*i*(-r13 + r31) + G3*r33*R31 + G4*r44*R41;
-						dr12_dt = -(gt*r12) + i*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1*r32);
-						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1*r11) - E2*r12 - delta1*r13 + E1*r33);
-						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E3*r12) - (delta1 - delta2 + delta3)*r14 + E1*r34);
-						//dr21_dt = -(gt*r21) + i*((delta1 - delta2)*r21 - E1*r23 + E2*r31 + E3*r41);
-						dr22_dt = gt/2. - gt*r22 + i*(-(E2*r23) - E3*r24 + E2*r32 + E3*r42) + G3*r33*R32 + G4*r44*R42;
-						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1*r21) - E2*r22 - delta2*r23 + E2*r33 + E3*r43);
-						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E3*r22) - delta3*r24 + E2*r34 + E3*r44);
-						//dr31_dt = -((G3 + 2*gt)*r31)/2. + i*(E1*r11 + E2*r21 + delta1*r31 - E1*r33);
-						//dr32_dt = -((G3 + 2*gt)*r32)/2. + i*(E1*r12 + E2*r22 + delta2*r32 - E2*r33 - E3*r34);
-						dr33_dt = i*(E1*r13 + E2*r23 - E1*r31 - E2*r32) - (G3 + gt)*r33;
-						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1*r14 + E2*r24 - E3*r32 + (delta2 - delta3)*r34);
-						//dr41_dt = -((G4 + 2*gt)*r41)/2. + i*(E3*r21 + (delta1 - delta2 + delta3)*r41 - E1*r43);
-						//dr42_dt = -((G4 + 2*gt)*r42)/2. + i*(E3*r22 + delta3*r42 - E2*r43 - E3*r44);
-						//dr43_dt = -((G3 + G4 + 2*gt)*r43)/2. + i*(E3*r23 - E1*r41 - E2*r42 + (-delta2 + delta3)*r43);
-						//dr44_dt = E3*i*(r24 - r42) - (G4 + gt)*r44;
+						dr11_dt = gt/2. - gt*r11 + G3*r33*R31 + G4*r44*R41 + i*(-(E1*r13) + E1c*r31);
+						dr12_dt = -(gt*r12) + i*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1c*r32);
+						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(delta1*r13) - E1c*r11 - E2c*r12 + E1c*r33);
+						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-((delta1 - delta2 + delta3)*r14) - E3c*r12 + E1c*r34);
+						dr22_dt = gt/2. - gt*r22 + G3*r33*R32 + G4*r44*R42 + i*(-(E2*r23) - E3*r24 + E2c*r32 + E3c*r42);
+						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(delta2*r23) - E1c*r21 - E2c*r22 + E2c*r33 + E3c*r43);
+						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(delta3*r24) - E3c*r22 + E2c*r34 + E3c*r44);
+						dr33_dt = -((G3 + gt)*r33) + i*(E1*r13 + E2*r23 - E1c*r31 - E2c*r32);
+						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*((delta2 - delta3)*r34 + E1*r14 + E2*r24 - E3c*r32);
+						dr44_dt = -((G4 + gt)*r44) + i*(E3*r24 - E3c*r42);
 					]]>
         </operator>
 				<operator kind="ex" constant="yes">