xmds2 model of Gena's 4-level system.

GenerateGenasSystem.nb: generate evolution
 equations
GenasSystemPlots: code for reading output data,
  converting complex field amplitudes to
  rotation angle/ellipticity, and plotting

1 files changed, 218 insertions, 0 deletions

diff --git a/xmds2/Genas_system/Genas_system.xmds b/xmds2/Genas_system/Genas_system.xmds
new file mode 100755
index 0000000..e855e4f
--- /dev/null
+++ b/xmds2/Genas_system/Genas_system.xmds
@@ -0,0 +1,218 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+	<name>Genas_system</name>
+	
+	<author>Eugeniy Mikhailov and Simon Rochester</author>
+	<description>
+		License GPL.
+
+		Solving 4 level atom in 0->1 configuration,
+		with field propagation along spatial axis Z
+		no Doppler broadening
+
+		We are solving 
+			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 - 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 c=3.e8;
+				const double lambda=794.7e-9; //wavelength in m
+				const double N=1e6*(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/16.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 [1/s]
+				const double gt=0.01*(2*M_PI*1e6);
+				// Natural linewidth  of upper state in [1/s]
+				const double g0=1.*(2*M_PI*1e6);
+
+				complex  Exc, Eyc; // Complex-conjugated Rabi frequency
+
+				complex  r21, r31, r41, r32, r42, r43, r44;   // density matrix elements
+			]]>
+		</globals>
+		<benchmark />
+		<arguments>
+			<!-- Real and imaginary parts of complex Rabi frequency in [1/s] -->
+			<argument name="ExReo" type="real" default_value="(3.+0.001)*(2*M_PI*1.e6)" />
+			<argument name="ExImo" type="real" default_value="0." />
+			<argument name="EyReo" type="real" default_value="0." />
+			<argument name="EyImo" type="real" default_value="(3.-0.001)*(2*M_PI*1.e6)" />
+			<!-- light detuning in [1/s] -->
+			<argument name="delta"  type="real" default_value="3.0*(2*M_PI*1e6)" />
+			<!--shift of upper-state M=0 sublevel-->
+			<argument name="delta0" type="real" default_value="1.*(2*M_PI*1e6)" />
+			<!--Static B-field Larmor frequency-->
+			<argument name="OL" type="real" default_value="4.*(2*M_PI*1e6)" />
+			<!--rf Rabi frequency-->
+			<argument name="Orf" type="real" default_value="0.1*(2*M_PI*1e6)" />
+			<!--rf frequency-->
+			<argument name="orf" type="real" default_value="4.*(2*M_PI*1e6)" />	
+		</arguments>
+		<bing />
+		<fftw plan="patient" />
+		<openmp />
+		<auto_vectorise />
+	</features>
+
+	<!-- 'z' and 't' to have dimensions [m] and [s]   -->
+	<geometry>
+		<propagation_dimension> z </propagation_dimension>
+		<transverse_dimensions>
+			<dimension name="t" lattice="1000" domain="(-1.5e-7, 2.5e-7)" />
+		</transverse_dimensions>
+	</geometry>
+
+	<!-- Rabi frequency --> 
+	<vector name="E_field" type="complex" initial_space="t">
+		<components>Ex Ey</components>
+		<initialisation>
+			<![CDATA[
+				// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
+				// as well as time
+				Ex=(ExReo+i*ExImo)*exp(-pow( ((t-0.0)/1e-7),2) );
+				Ey=(EyReo+i*EyImo)*exp(-pow( ((t-0.0)/1e-7),2) );
+			]]>
+		</initialisation>
+	</vector>
+
+	<vector name="density_matrix" type="complex" initial_space="t">
+		<components>r11 r22 r33 r12 r13 r14 r23 r24 r34 r44</components>
+		<initialisation>
+			<![CDATA[
+				r11 = 1; r22 = 0; r33 = 0; r44 = 0;
+				r12 = 0; r13 = 0; r14 = 0;
+				r23 = 0; r24 = 0;
+				r34 = 0;
+			]]>
+		</initialisation>
+	</vector>
+
+	<vector name="rfField" type="real">
+		<components> rf </components>
+		<initialisation>
+		  <![CDATA[
+			rf = Orf*sin(orf*t);
+		  ]]>
+		</initialisation>
+	</vector>	
+	
+	<sequence>
+		<integrate algorithm="ARK45" tolerance="0.05e-7" interval="1e1">
+			<samples>200 200</samples>
+			<operators>
+				<operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
+					<integration_vectors>density_matrix</integration_vectors>
+					<dependencies>E_field rfField</dependencies>
+					<boundary_condition kind="left">
+						<![CDATA[
+							r11 = 1; r22 = 0; r33 = 0; r44 = 0;
+							r12 = 0; r13 = 0; r14 = 0;
+							r23 = 0; r24 = 0;
+							r34 = 0;
+						]]>
+					</boundary_condition>
+					<![CDATA[
+						Exc = conj(Ex);
+						Eyc = conj(Ey);
+
+						r21=conj(r12);
+						r31=conj(r13);
+						r41=conj(r14);
+						r32=conj(r23);
+						r42=conj(r24);
+						r43=conj(r34);
+
+						// Equations of motions according to Simon's mathematica code
+						dr11_dt = gt - gt*r11 + (-Ey - Ex*i)*r12 + i*(Ex + Ey*i)*r14 + i*(Exc + Eyc*i)*r21 + (-Eyc - Exc*i)*r41 + g0*(r22 + r33 + r44);
+						dr12_dt = (2*(Eyc - Exc*i)*r11 - i*((2*delta - (g0 + 2*gt)*i - 2*OL - 2*rf)*r12 - 2*(Exc + Eyc*i)*r22 + 2*(Exc - Eyc*i)*r42))/2.;
+						dr13_dt = (-g0/2. - gt - delta*i + delta0*i)*r13 + i*(Exc + Eyc*i)*r23 + (-Eyc - Exc*i)*r43;
+						dr14_dt = -(i*((2*delta - (g0 + 2*gt)*i + 2*OL + 2*rf)*r14 - 2*(Exc + Eyc*i)*r24 - 2*(Exc - Eyc*i)*(r11 - r44)))/2.;
+						dr22_dt = (Ey + Ex*i)*r12 + (Eyc - Exc*i)*r21 - (g0 + gt)*r22;
+						dr23_dt = (Ey + Ex*i)*r13 + i*(delta0 + i*(g0 + gt + i*OL) - rf)*r23;
+						dr24_dt = (Ey + Ex*i)*r14 + (Eyc + Exc*i)*r21 - (g0 + gt + 2*i*OL + 2*i*rf)*r24;
+						dr33_dt = -((g0 + gt)*r33);
+						dr34_dt = (Eyc + Exc*i)*r31 - i*(delta0 - (g0 + gt)*i + OL + rf)*r34;
+						dr44_dt = (Ey - Ex*i)*r14 + (Eyc + Exc*i)*r41 - (g0 + gt)*r44;
+					]]>
+				</operator>
+				<operator kind="ex" constant="yes">
+					<operator_names>Lt</operator_names>
+					<![CDATA[
+						Lt = i*1./c*kt;
+					]]>
+				</operator>
+				<integration_vectors>E_field</integration_vectors>
+				<dependencies>density_matrix</dependencies>
+				<![CDATA[
+					dEx_dz = i*eta*conj(r12-r14) - Lt[Ex] ;
+					dEy_dz = -eta*conj(r12+r14) - Lt[Ey] ;
+				]]>
+			</operators>
+		</integrate>
+	</sequence>
+
+	<!-- The output to generate -->
+	<output format="binary" filename="Genas_system.xsil">
+		<group>
+			<sampling basis="t(100)" initial_sample="yes">
+				<dependencies>E_field</dependencies>
+				<moments>Ex_re_out Ex_im_out Ey_re_out Ey_im_out</moments>
+				<![CDATA[
+					Ex_re_out = Ex.Re();
+					Ex_im_out = Ex.Im();
+					Ey_re_out = Ey.Re();
+					Ey_im_out = Ey.Im();
+				]]>
+			</sampling>
+		</group>
+
+		<group>
+			<sampling basis="t(100)" initial_sample="yes">
+				<dependencies>density_matrix</dependencies>
+				<moments>
+					r11_out r22_out r33_out r44_out 
+					r12_re_out r12_im_out r13_re_out r13_im_out r14_re_out r14_im_out
+					                      r23_re_out r23_im_out r24_re_out r24_im_out 
+					                                            r34_re_out r34_im_out
+				</moments>
+				<![CDATA[
+					// populations output 
+					r11_out = r11.Re();
+					r22_out = r22.Re();
+					r33_out = r33.Re();
+					r44_out = r44.Re();
+					// coherences output 
+					r12_re_out = r12.Re();
+					r12_im_out = r12.Im();
+					r13_re_out = r13.Re();
+					r13_im_out = r13.Im();
+					r14_re_out = r14.Re();
+					r14_im_out = r14.Im();
+					r23_re_out = r23.Re();
+					r23_im_out = r23.Im();
+					r24_re_out = r24.Re();
+					r24_im_out = r24.Im();
+					r34_re_out = r34.Re();
+					r34_im_out = r34.Im();
+				]]>
+			</sampling>
+		</group>
+	</output>
+</simulation>
+	
+<!-- 
+vim: ts=2 sw=2 foldmethod=indent: 
+-->