matches the Simon;s mathematica code

1 files changed, 36 insertions, 32 deletions

diff --git a/xmds2/Nlevels_no_dopler_no_z.example/Nlevels_no_dopler_no_z_from_mathematica.xmds b/xmds2/Nlevels_no_dopler_no_z.example/Nlevels_no_dopler_no_z_from_mathematica.xmds
index 04c7306..3f306c5 100644
--- a/xmds2/Nlevels_no_dopler_no_z.example/Nlevels_no_dopler_no_z_from_mathematica.xmds
+++ b/xmds2/Nlevels_no_dopler_no_z.example/Nlevels_no_dopler_no_z_from_mathematica.xmds
@@ -51,8 +51,8 @@
 			// repopulation rate (atoms flying in/out the laser beam)  in MHz
 			const double gt=0.01/2;
 			// Natural linewidth  of j's level in MHz
-			const double G3=2.7;
-			const double G4=3.0;
+			const double G3=2*2.7;
+			const double G4=2*3.0;
 
 			// branching ratios
 			const double R41=0.5, R42=0.5;
@@ -68,7 +68,7 @@
 
 			complex  E1c, E2c, Ecc; // Complex conjugated Rabi frequencies
 
-			//complex  r21, r31, r41, r32, r42, r43;
+			complex  r21, r31, r41, r32, r42, r43, r44;
 
 			]]>
 		</globals>
@@ -82,7 +82,9 @@
 	</geometry>
 
 	<vector name="density_matrix" type="complex" dimensions="d1">
-		<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   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>
 		<!--
 				 note one of the level population is redundant since
 				 r11+r22+r33+r44=1
@@ -92,14 +94,14 @@
 			<![CDATA[
 			// Note: 
 			// convergence is really slow if all populations concentrated at the bottom level |1>
-			// this is because if r11=1, everything else is 0 and then every small incriment 
+			// this is because if r11=1, everything else is 0 and then every small increment 
 			// seems to be huge and adaptive solver makes smaller and smaller steps.
-			// As quicj and dirty fix I reshufle initial poulation  
+			// As quick and dirty fix I reshuffle initial population  
 			// so some of the population sits at the  second ground level |2>
 			// TODO: Fix above. Make the equation of motion for r11 
 			//       and express other level, let's say r44
 			//       through population normalization
-			r22 = 0.001; r33 = 0; r44 = 0;
+			r11 = 1; r22 = 0; r33 = 0;
 			r12 = 0; r13 = 0; r14 = 0;
 			r23 = 0; r24 = 0;
 			r34 = 0;
@@ -108,7 +110,7 @@
 	</vector>
 
 	<sequence>
-		<integrate algorithm="ARK45" interval="600" tolerance="1e-3">
+		<integrate algorithm="ARK89" interval="3" tolerance="1e-3">
 			<samples>128</samples>
 			<operators>
 				<integration_vectors>density_matrix</integration_vectors>
@@ -117,30 +119,32 @@
 				E2c = conj(E2);
 				Ecc = conj(Ec);
 
-				// r21=conj(r12);
-				// r31=conj(r13);
-				// r41=conj(r14);
-				// r32=conj(r23);
-				// r42=conj(r24);
-				// r43=conj(r34);
-
-				// Equations of motions see page 1928 of the Frank's JMO paper
-				dr11_dt = gt - 2*gt*r11 + E1*i*(-r13 + r31) + 2*G3*r33*R31 + 2*G4*r44*R41;
-				dr12_dt = -2*gt*r12 + i*((-d1 + d2)*r12 - E2*r13 - Ec*r14 + E1*r32);
-				dr13_dt = -((G3 + 2*gt)*r13) + i*(-(E1*r11) - E2*r12 - d1*r13 + E1*r33);
-				dr14_dt = -((G4 + 2*gt)*r14) + i*(-(Ec*r12) - (d1 - d2 + d3)*r14 + E1*r34);
-				dr21_dt = -2*gt*r21 + i*((d1 - d2)*r21 - E1*r23 + E2*r31 + Ec*r41);
-				dr22_dt = gt - 2*gt*r22 + i*(-(E2*r23) - Ec*r24 + E2*r32 + Ec*r42) + 2*G3*r33*R32 + 2*G4*r44*R42;
-				dr23_dt = -((G3 + 2*gt)*r23) + i*(-(E1*r21) - E2*r22 - d2*r23 + E2*r33 + Ec*r43);
-				dr24_dt = -((G4 + 2*gt)*r24) + i*(-(Ec*r22) - d3*r24 + E2*r34 + Ec*r44);
-				dr31_dt = -((G3 + 2*gt)*r31) + i*(E1*r11 + E2*r21 + d1*r31 - E1*r33);
-				dr32_dt = -((G3 + 2*gt)*r32) + i*(E1*r12 + E2*r22 + d2*r32 - E2*r33 - Ec*r34);
-				dr33_dt = i*(E1*r13 + E2*r23 - E1*r31 - E2*r32) - 2*(G3 + gt)*r33;
-				dr34_dt = -((G3 + G4 + 2*gt)*r34) + i*(E1*r14 + E2*r24 - Ec*r32 + (d2 - d3)*r34);
-				dr41_dt = -((G4 + 2*gt)*r41) + i*(Ec*r21 + (d1 - d2 + d3)*r41 - E1*r43);
-				dr42_dt = -((G4 + 2*gt)*r42) + i*(Ec*r22 + d3*r42 - E2*r43 - Ec*r44);
-				dr43_dt = -((G3 + G4 + 2*gt)*r43) + i*(Ec*r23 - E1*r41 - E2*r42 + (-d2 + d3)*r43);
-				dr44_dt = Ec*i*(r24 - r42) - 2*(G4 + gt)*r44;
+				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 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*((-d1 + d2)*r12 - E2*r13 - Ec*r14 + E1*r32);
+				dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1*r11) - E2*r12 - d1*r13 + E1*r33);
+				dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(Ec*r12) - (d1 - d2 + d3)*r14 + E1*r34);
+				//dr21_dt = -(gt*r21) + i*((d1 - d2)*r21 - E1*r23 + E2*r31 + Ec*r41);
+				dr22_dt = gt/2. - gt*r22 + i*(-(E2*r23) - Ec*r24 + E2*r32 + Ec*r42) + G3*r33*R32 + G4*r44*R42;
+				dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1*r21) - E2*r22 - d2*r23 + E2*r33 + Ec*r43);
+				dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(Ec*r22) - d3*r24 + E2*r34 + Ec*r44);
+				//dr31_dt = -((G3 + 2*gt)*r31)/2. + i*(E1*r11 + E2*r21 + d1*r31 - E1*r33);
+				//dr32_dt = -((G3 + 2*gt)*r32)/2. + i*(E1*r12 + E2*r22 + d2*r32 - E2*r33 - Ec*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 - Ec*r32 + (d2 - d3)*r34);
+				//dr41_dt = -((G4 + 2*gt)*r41)/2. + i*(Ec*r21 + (d1 - d2 + d3)*r41 - E1*r43);
+				//dr42_dt = -((G4 + 2*gt)*r42)/2. + i*(Ec*r22 + d3*r42 - E2*r43 - Ec*r44);
+				//dr43_dt = -((G3 + G4 + 2*gt)*r43)/2. + i*(Ec*r23 - E1*r41 - E2*r42 + (-d2 + d3)*r43);
+				//dr44_dt = Ec*i*(r24 - r42) - (G4 + gt)*r44;
+
 				]]>
 			</operators>
 		</integrate>