First draft of Doppler averaging

1 files changed, 137 insertions, 71 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 ae2fe75..7d33920 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
@@ -45,11 +45,23 @@
 			<![CDATA[
 				const double pi = M_PI; 
 				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
 				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)]
+				
+				//  ---------  Atom and cell properties -------------------------
+				// range of Maxwell distribution atomic velocities
+				const double Temperature=1000.01; // cell temperature in [K] ! make sure it is not equal to zero!
+				const double mass = (86.909180527 * 1.660538921e-27); // atom mass in [kg] 
+				// above mass expression is written as (expression is isotopic_mass * atomic_mass_unit)
+
+				// Average sqrt(v^2) in Maxwell distribution for one dimension
+				const double v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell distribution velocities range to take in account in [m/s]
+				const double V_maxwell_min = -10, V_maxwell_max = -V_maxwell_min; // TODO: express via mass and Temperature
 
 				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
 				const double gt=0.01 *(2*M_PI*1e6);
@@ -62,12 +74,16 @@
 				const double R31=0.5, R32=0.5;
 
 
-				complex  E1c, E2c, E3c, E4c; // Complex conjugated Rabi frequencies
+				complex  E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies
 
 				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
 
+				// inner use variables 
+				double probability_v; // will be used as p(v) in Maxwell distribution
+
 			]]>
 		</globals>
+		<validation kind="run-time"/> <!--allows to put ranges as wariables-->
 		<benchmark />
     <arguments>
 			<!-- Rabi frequency divided by 2 in [1/s] -->
@@ -91,7 +107,7 @@
 		<propagation_dimension> z </propagation_dimension>
 		<transverse_dimensions>
 			<dimension name="t"   lattice="1000"   domain="(-14.0e-7, 4.0e-7)" />
-			<dimension name="v"   lattice="10"   domain="(-1000, 1000)" />
+			<dimension name="v"   lattice="10"   domain="(V_maxwell_min, V_maxwell_max)" />
 		</transverse_dimensions>
 	</geometry>
 
@@ -110,17 +126,44 @@
 		</initialisation>
 	</vector>
 
-	<!--
-	<computed_vector name="E_field_avgd" dimensions="z t" type="complex">
+	<!--Maxwell distribution probability p(v)-->
+	<computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real">
+		<components>probability_v</components>
+		<evaluation>
+			<![CDATA[
+			probability_v=1.0/(v_thermal_averaged*sqrt(2*M_PI)) * exp( - mod2(v/v_thermal_averaged)/2.0 ); 
+			]]>
+		</evaluation>
+	</computed_vector>
+
+	<!--Maxwell distribution norm sum(p(v))
+			 Needed since we sum over the grid instead of true integral,
+			 we also have finite cut off velocities-->
+	<computed_vector name="Maxwell_distribution_probabilities_norm" dimensions="" type="real">
+		<components>probability_v_norm</components>
+		<evaluation>
+			<dependencies basis="v">Maxwell_distribution_probabilities</dependencies>
+			<![CDATA[
+			probability_v_norm=probability_v;
+			]]>
+		</evaluation>
+	</computed_vector>
+
+
+	<!-- Averaged across Maxwell distribution fields amplitudes -->
+	<computed_vector name="E_field_avgd" dimensions="t" type="complex">
 		<components>E1a E2a E3a E4a</components>
 		<evaluation>
-			<dependencies basis="v">E_field</dependencies>
+			<dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
 			<![CDATA[
-			// Calculate the current normalisation of the wave function.
+			double prob_v_normalized=probability_v/probability_v_norm;
+			E1a=E1*prob_v_normalized;
+			E2a=E2*prob_v_normalized;
+			E3a=E3*prob_v_normalized;
+			E4a=E4*prob_v_normalized;
 			]]>
 		</evaluation>
 	</computed_vector>
-	-->
 
 
 	<vector name="density_matrix" type="complex" initial_space="t">
@@ -155,55 +198,55 @@
 			<!-- 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</dependencies>
+				<dependencies>E_field_avgd</dependencies>
 				<![CDATA[
-						E1c = conj(E1);
-						E2c = conj(E2);
-						E3c = conj(E3);
-						E4c = conj(E4);
+						E1ac = conj(E1a);
+						E2ac = conj(E2a);
+						E3ac = conj(E3a);
+						E4ac = conj(E4a);
 
 				// IMPORTANT: assumes no detunings
-				r11 = (gt*(4*mod2(E1) + (G3 + gt)*(G3 + 2*gt))*((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt)) + 4*mod2(E3)*G4*(R41 -
-				R42)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r13 = -((E1c*gt*(G3 + gt)*i*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt))
-				+ 4*mod2(E3)*G4*(R41 -
-				R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r22 = (gt*(4*mod2(E3) + (G4 + gt)*(G4 + 2*gt))*((G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt)) + 4*mod2(E1)*G3*(-R31 +
-				R32)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r24 = -((E3c*gt*(G4 + gt)*i*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt))
-				+ 4*mod2(E1)*G3*(-R31 +
-				R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
-
-				r33 = (2*mod2(E1)*gt*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt)) +
-				4*mod2(E3)*G4*(R41 -
-				R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42));
-
-				r44 = (2*mod2(E3)*gt*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt)) +
-				4*mod2(E1)*G3*(-R31 +
-				R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
-				2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
-				gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
-				2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42));
+				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>
@@ -215,16 +258,26 @@
 		<!--<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">
+		<integrate algorithm="RK4" steps="100" tolerance="1e-5" interval="1.5e-3">
 		<!--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-->
 
 		<!--<integrate algorithm="SIC" interval="4e-2" steps="200">-->
-			<samples>200 200</samples>
+
+			<filters where="step end">
+				<filter>
+					<dependencies>E_field_avgd</dependencies>
+					<![CDATA[
+					//printf("z=%g, t=%g\n", z, t);
+					//printf("E_field_avgd filter in integrator\n", z, t); 
+					]]>
+				</filter>
+			</filters>
+			<samples>100 100 100</samples>
 			<operators>
         <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
 					<integration_vectors>density_matrix</integration_vectors>
-          <dependencies>E_field</dependencies>
+          <dependencies>E_field_avgd</dependencies>
           <boundary_condition kind="left">
 						<!--This set boundary condition at all 'z'  and left border of 't' (i.e. min(t))-->
 						<!--
@@ -238,10 +291,10 @@
 						-->
           </boundary_condition>
 					<![CDATA[
-						E1c = conj(E1);
-						E2c = conj(E2);
-						E3c = conj(E3);
-						E4c = conj(E4);
+						E1ac = conj(E1a);
+						E2ac = conj(E2a);
+						E3ac = conj(E3a);
+						E4ac = conj(E4a);
 
 						r21=conj(r12);
 						r31=conj(r13);
@@ -252,17 +305,17 @@
 						//r44=1- r11 - r22 - r33;
 
 						// 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 + 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 + 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 + Kvec*v) - (delta3 + Kvec*v))*r34);
-						dr44_dt = i*(E4*r14 + E3*r24 - E4c*r41 - E3c*r42) - (G4 + gt)*r44;
+						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;
 					]]>
         </operator>
 				<operator kind="ex" constant="yes">
@@ -331,6 +384,19 @@
 				]]>
 			</sampling>
 		</group>
+
+		<group>
+      <sampling basis="t(100) " initial_sample="yes">
+				<dependencies>E_field_avgd</dependencies>
+				<moments>I1_out_avgd I2_out_avgd I3_out_avgd I4_out_avgd</moments>
+				<![CDATA[
+				I1_out_avgd = mod2(E1a);
+				I2_out_avgd = mod2(E2a);
+				I3_out_avgd = mod2(E3a);
+				I4_out_avgd = mod2(E4a);
+				]]>
+			</sampling>
+		</group>
 	</output>
 
 </simulation>