Eta parameter fixes, Rabi measured in [1/s] and z in [m]

Very tedious check of units, now everything in SI units
which is a bit painful for Rabi frequencies.

1 files changed, 30 insertions, 39 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 718cf76..a210522 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
@@ -29,62 +29,52 @@
 		*                    ------- |1>
 		*
 
-		We moved to dimensionless units
-		t -> t*DecayRateNormalization    ,time
-		z -> z*DecayRateNormalization/c  , distance
-		rabi_frequency -> rabi_frequency/DecayRateNormalization 
-		eta -> eta*c/DecayRateNormalization^2      , coupling constant
-		gij -> gij/DecayRateNormalization
-		Wij -> Wij/DecayRateNormalization
-
-		where g is 1MHz rate
 
 		We are solving 
-			dE/dz+(1/c)*dE/dt=i*eta*rho_ij,  where j llevel is higher then i
+			dE/dz+(1/c)*dE/dt=i*eta*rho_ij,  where j level is higher then i.
+			Note that E is actually a Rabi frequency of electromagnetic field not the EM field
 		in xmds terms it looks like
-			dE_dz = i*eta*rhoij - L[E], here we moved t dependence to Fourier space
+			dE_dz = i*eta*rhoij - 1/c*L[E], here we moved t dependence to Fourier space
+
+		VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to
+		normalize it to something else look drho/dt equation.
+		No need to renormalizes eta as long as its express through i
+		the upper level decay rate in the same units as Rabi frequency.
 	</description>
 
 	<features>
 		<globals>
 			<![CDATA[
 				const double pi = M_PI; 
-				const double DecayRateNormalization=10e6; // measured in s^-1
-				const double c=3e8;
+				const double c=3.e8;
 				const double lambda=794.7e-9; //wavelength in m
-				const double N=1e9; //number of particles per cubic cm
-				const double eta = 3*lambda*lambda*(N/1e-6)/8.0/pi*c/DecayRateNormalization;  //dimensionless coupling constant
+				const double N=1e10*(1e6);    //number of particles per cubic m i.e. density
+				const double Gamma=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/8.0/M_PI;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
 
-				// 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*2.7;
-				const double G4=2*3.0;
+				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
+				const double gt=0.01/2 *(2*M_PI*1e6);
+				// Natural linewidth  of j's level in [1/s]
+				const double G3=2*2.7 *(2*M_PI*1e6);
+				const double G4=2*3.0 *(2*M_PI*1e6);
 
-				// branching ratios
+				// 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;
 
-				//const double delta1=0;         // E2 detuning in MHz
-				//const double delta2=0;         // E2 detuning in MHz
-				//const double delta3=0;         // E3 detuning in MHz
-
-				//const complex E1o=0.005/2;     // Rabi frequency in MHz
-				//const complex E2o=5.0/2;         // Rabi frequency in MHz
-				//const complex E3o=6.0/2;         // Rabi frequency in MHz
 
 				complex  E1c, E2c, E3c; // Complex conjugated Rabi frequencies
 
-				complex  r21, r31, r41, r32, r42, r43, r44;
+				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
 			]]>
 		</globals>
 		<benchmark />
     <arguments>
-			<!--Rabi frequency divided by 2-->
-      <argument name="E1o" type="real" default_value="0.0025" />
-      <argument name="E2o" type="real" default_value="2.5" />
-      <argument name="E3o" type="real" default_value="3.0" />
-			<!--Fields detuning in MHz-->
+			<!-- Rabi frequency divided by 2 in [1/s] -->
+      <argument name="E1o" type="real" default_value="0.0025*(2*M_PI*1e6)" />
+      <argument name="E2o" type="real" default_value="6*(2*M_PI*1e6)" />
+      <argument name="E3o" type="real" default_value="0.0*(2*M_PI*1e6)" />
+			<!-- Fields detuning in [1/s] -->
       <argument name="delta1"  type="real" default_value="0.0" />
       <argument name="delta2"  type="real" default_value="0.0" />
       <argument name="delta3"  type="real" default_value="0.0" />
@@ -95,21 +85,22 @@
 		<auto_vectorise />
 	</features>
 
-	<!-- Let's convert 'z' and 't' to dimentionless units over here -->
+	<!-- 'z' and 't' to have dimensions [m] and [s]   -->
 	<geometry>
 		<propagation_dimension> z </propagation_dimension>
 		<transverse_dimensions>
-			<dimension name="t"   lattice="1000"   domain="(-2.0, 4.0)" />
+			<dimension name="t"   lattice="1000"   domain="(-2.0e-6, 4.0e-6)" />
 		</transverse_dimensions>
 	</geometry>
 
+	<!-- Rabi frequency --> 
 	<vector name="E_field" type="complex" initial_space="t">
 		<components>E1 E2 E3</components>
 		<initialisation>
 			<![CDATA[
 			// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
 			// as well as time
-			E1=E1o*exp(-pow( ((t-0.0)/1),2) );
+			E1=E1o*exp(-pow( ((t-0.0)/1e-6),2) );
 			E2=E2o;
 			E3=E3o;
 			]]>
@@ -146,7 +137,7 @@
 
 	<sequence>
 		<!--For this set of conditions ARK45 is faster than ARK89-->
-		<integrate algorithm="ARK45" tolerance="1e-2" interval="4e-2">
+		<integrate algorithm="ARK45" tolerance="1e-5" interval="7e-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-->
 
@@ -193,7 +184,7 @@
 				<operator kind="ex" constant="yes">
 					<operator_names>Lt</operator_names>
 					<![CDATA[
-						Lt = i*kt;
+					Lt = i*1./c*kt;
 					]]>
 				</operator>
         <integration_vectors>E_field</integration_vectors>