Added velocity dimension

need to do the doppler averaging

1 files changed, 29 insertions, 11 deletions

diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds
index 8faa04d..ae2fe75 100644
--- a/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds
+++ b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds
@@ -46,6 +46,7 @@
 				const double pi = M_PI; 
 				const double c=3.e8;
 				const double lambda=794.7e-9; //wavelength in m
+				const double Kvec = 2*M_PI/lambda;  // k-vector
 				const double N=1e09*(1e6);    //number of particles per cubic m i.e. density
 				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
 				const double eta = 3*lambda*lambda*N*Gamma_super/8.0/M_PI;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
@@ -64,6 +65,7 @@
 				complex  E1c, E2c, E3c, E4c; // Complex conjugated Rabi frequencies
 
 				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
+
 			]]>
 		</globals>
 		<benchmark />
@@ -89,7 +91,7 @@
 		<propagation_dimension> z </propagation_dimension>
 		<transverse_dimensions>
 			<dimension name="t"   lattice="1000"   domain="(-14.0e-7, 4.0e-7)" />
-			<dimension name="v"   lattice="1"   domain="(-1000, 1000)" />
+			<dimension name="v"   lattice="10"   domain="(-1000, 1000)" />
 		</transverse_dimensions>
 	</geometry>
 
@@ -108,6 +110,19 @@
 		</initialisation>
 	</vector>
 
+	<!--
+	<computed_vector name="E_field_avgd" dimensions="z t" type="complex">
+		<components>E1a E2a E3a E4a</components>
+		<evaluation>
+			<dependencies basis="v">E_field</dependencies>
+			<![CDATA[
+			// Calculate the current normalisation of the wave function.
+			]]>
+		</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>-->
@@ -196,7 +211,11 @@
 
 	<sequence>
 		<!--For this set of conditions ARK45 is faster than ARK89-->
-		<integrate algorithm="ARK45" tolerance="1e-5" interval="1.5e-2">
+		<!--ARK45 is good for small detuning when all frequency like term are close to zero-->
+		<!--<integrate algorithm="ARK45" tolerance="1e-5" interval="1.5e-2">--> 
+		<!--RK4 is good for large detunings when frequency like term are big, it does not try to be too smart about adaptive step which ARK seems to make too small-->
+		<!--When ARK45 works it about 3 times faster then RK4 with 1000 steps-->
+		<integrate algorithm="RK4" steps="1000" tolerance="1e-5" interval="1.5e-2">
 		<!--SIC algorithm seems to be much slower and needs fine 'z'  step tuning and much finer time grid-->
 		<!--For example I had to quadruple the time grid from 1000 to 4000 when increased z distance from 0.02 to 0.04-->
 
@@ -234,15 +253,15 @@
 
 						// Equations of motions according to Nate's mathematica code
 						dr11_dt = gt/2 - gt*r11 + i*(-(E1*r13) - E4*r14 + E1c*r31 + E4c*r41) + G3*r33*R31 + G4*r44*R41;
-						dr12_dt = -(gt*r12) + i*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1c*r32 + E4c*r42);
-						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1c*r11) - E2c*r12 - delta1*r13 + E1c*r33 + E4c*r43);
-						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4c*r11) - E3c*r12 - (delta1 - delta2 + delta3)*r14 + E1c*r34 + E4c*r44);
+						dr12_dt = -(gt*r12) + i*((-(delta1 + Kvec*v) + (delta2 + Kvec*v))*r12 - E2*r13 - E3*r14 + E1c*r32 + E4c*r42);
+						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1c*r11) - E2c*r12 - (delta1 + Kvec*v)*r13 + E1c*r33 + E4c*r43);
+						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4c*r11) - E3c*r12 - ((delta1 + Kvec*v) - (delta2 + Kvec*v) + (delta3 + Kvec*v))*r14 + E1c*r34 + E4c*r44);
 						dr22_dt = gt/2 - gt*r22 + i*(-(E2*r23) - E3*r24 + E2c*r32 + E3c*r42) + G3*r33*R32 + G4*r44*R42;
-						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1c*r21) - E2c*r22 - delta2*r23 + E2c*r33 + E3c*r43);
-						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4c*r21) - E3c*r22 - delta3*r24 + E2c*r34 + E3c*r44);
+						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1c*r21) - E2c*r22 - (delta2 + Kvec*v)*r23 + E2c*r33 + E3c*r43);
+						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4c*r21) - E3c*r22 - (delta3 + Kvec*v)*r24 + E2c*r34 + E3c*r44);
 						dr33_dt = i*(E1*r13 + E2*r23 - E1c*r31 - E2c*r32) - (G3 + gt)*r33;
 
-						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1*r14 + E2*r24 - E4c*r31 - E3c*r32 + (delta2 - delta3)*r34);
+						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*(E1*r14 + E2*r24 - E4c*r31 - E3c*r32 + ((delta2 + Kvec*v) - (delta3 + Kvec*v))*r34);
 						dr44_dt = i*(E4*r14 + E3*r24 - E4c*r41 - E3c*r42) - (G4 + gt)*r44;
 					]]>
         </operator>
@@ -266,11 +285,10 @@
 
 
 
-
 	<!-- The output to generate -->
 	<output format="binary" filename="Nlevels_with_doppler_with_z_4wm.xsil">
 		<group>
-      <sampling basis="t(1000)" initial_sample="yes">
+      <sampling basis="t(100) v(10)" initial_sample="yes">
 				<dependencies>E_field</dependencies>
 				<moments>I1_out I2_out I3_out I4_out</moments>
 				<![CDATA[
@@ -283,7 +301,7 @@
 		</group>
 
 		<group>
-      <sampling basis="t(100)" initial_sample="yes">
+      <sampling basis="t(100) v(10)" initial_sample="yes">
 				<dependencies>density_matrix</dependencies>
 				<moments>
 					r11_out r22_out r33_out r44_out