Nlevel example moved to more appropriate place

5 files changed, 440 insertions, 0 deletions

diff --git a/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Makefile b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Makefile
new file mode 100644
index 0000000..a89800a
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Makefile
@@ -0,0 +1,37 @@
+###  -*- make -*-
+### This file is part of the Debian xmds package
+### Copyright (C) 2006 Rafael Laboissiere
+### This file is relased under the GNU General Public License
+### NO WARRANTIES!
+
+### This makefile can be used to build and run the XMDS examples
+
+XMDS_FILES = $(shell ls *.xmds)
+RUN_FILES = $(patsubst %.xmds,%.run,$(XMDS_FILES))
+CC_FILES = $(patsubst %.xmds,%.cc,$(XMDS_FILES))
+XSIL_FILES = $(patsubst %.xmds,%.xsil,$(XMDS_FILES))
+M_FILES = $(patsubst %.xmds,%.m,$(XMDS_FILES))
+
+XMDS = xmds2
+XSIL2GRAPHICS = xsil2graphics
+
+all: $(M_FILES) 
+
+%.run: %.xmds
+	$(XMDS) $<
+	mv $(patsubst %.xmds,%,$<) $@
+
+%.xsil: %.run
+	./$<
+
+%.m: %.xsil
+	$(XSIL2GRAPHICS) $<
+
+plot: $(M_FILES)
+	octave pp.m
+
+clean:
+	rm -f $(CC_FILES) $(RUN_FILES) $(M_FILES) $(XSIL_FILES) *.wisdom.fftw3 *.dat octave-core *.wisdom *.pdf
+
+.PRECIOUS: %.run %.xsil %.m
+.PHONY: all clean
diff --git a/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Nlevels_no_dopler_no_z.xmds b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Nlevels_no_dopler_no_z.xmds
new file mode 100644
index 0000000..91b1711
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Nlevels_no_dopler_no_z.xmds
@@ -0,0 +1,179 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+	<name>Nlevels_no_dopler_no_z</name>
+
+	<author>Eugeniy Mikhailov</author>
+	<description>
+		License GPL.
+
+		Solving 4 level atom in N-field configuration,
+		no field propagation along spatial axis included
+		no Doppler broadening
+
+		For master equations look "Four-level 'N-scheme' in bare and quasi-dressed states pictures"
+		by T. Abi-Salloum, S. Meiselman, J.P. Davis and F.A. Narducci
+		Journal of Modern Optics, 56: 18, 1926 -- 1932, (2009).
+
+		Present calculation matches Fig.3 (b) from the above paper. 
+		Note that I need to double all Rabi frequencies to match the figure.
+
+		*  -------- |4>
+		*    \        
+		*     \ Ec    -------- |3>  
+		*      \       /    \
+		*       \  E2 /      \ 
+		*        \   /        \ E1
+		*       ------- |2>    \
+		*                       \
+		*                    ------- |1>
+		*
+
+		We moved to dimensionless units
+		t -> t*g    ,time
+		z -> z*g/c  , distance
+		rabi_frequency -> rabi_frequency/g 
+		eta -> eta*c/g^2      , coupling constant
+		gij -> gij/g
+		Wij -> Wij/g
+
+		where g is 1MHz rate
+	</description>
+
+	<features>
+		<benchmark />
+		<bing />
+		<fftw plan="patient" />
+		<openmp />
+		<auto_vectorise />
+		<globals>
+			<![CDATA[
+			// Wij decay rate of levels i-->j
+			const double W12=0.01;     // population decay rate in MHz
+			const double W13=0.0;      // population decay rate in MHz
+			const double W14=0.0;      // population decay rate in MHz
+			const double W21=W12;      // population decay rate in MHz
+			const double W23=0.0;      // population decay rate in MHz
+			const double W24=0.0;      // population decay rate in MHz
+			const double W31=2.7;      // population decay rate in MHz
+			const double W32=W31;      // population decay rate in MHz
+			const double W34=0.0;      // population decay rate in MHz
+			const double W41=3.0;      // population decay rate in MHz
+			const double W42=W41;      // population decay rate in MHz
+			const double W43=0.0;      // population decay rate in MHz
+
+			const double g12=(W12+W13+W14+W21+W23+W34)/2;  // coherence decay rate in MHz
+			const double g13=(W12+W13+W14+W31+W32+W34)/2;  // coherence decay rate in MHz
+			const double g14=(W12+W13+W14+W43+W42+W41)/2;  // coherence decay rate in MHz
+			const double g23=(W21+W23+W24+W31+W32+W34)/2;  // coherence decay rate in MHz
+			const double g24=(W21+W23+W24+W41+W42+W43)/2;  // coherence decay rate in MHz
+			const double g34=(W31+W32+W34+W41+W42+W43)/2;  // coherence decay rate in MHz
+
+			// const double d1=0;         // E2 detuning in MHz
+			const double d2=0;         // E2 detuning in MHz
+			const double dc=0;         // Ec detuning in MHz
+
+			const complex E1=0.005;     // Rabi frequency in MHz
+			const complex E2=5;         // Rabi frequency in MHz
+			const complex Ec=6;         // Rabi frequency in MHz
+
+			complex  E1c, E2c, Ecc; // Complex conjugated Rabi frequencies
+
+			complex r11, r21, r31, r41, r32, r42, r43;
+
+			]]>
+		</globals>
+	</features>
+
+	<geometry>
+		<propagation_dimension> t </propagation_dimension>
+		<transverse_dimensions>
+			<dimension name="d1" lattice="1280"  domain="(-60, 60)" />
+		</transverse_dimensions>
+	</geometry>
+
+	<vector name="density_matrix" type="complex" dimensions="d1">
+		<components>r22 r33 r44 r12 r13 r14 r23 r24 r34</components>
+		<!--
+				 note one of the level population is redundant since
+				 r11+r22+r33+r44=1
+				 so r11 is missing
+		-->
+		<initialisation>
+			<![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 increment 
+			// seems to be huge and adaptive solver makes smaller and smaller steps.
+			// 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;
+			r12 = 0; r13 = 0; r14 = 0;
+			r23 = 0; r24 = 0;
+			r34 = 0;
+			]]>
+		</initialisation>
+	</vector>
+
+	<sequence>
+		<integrate algorithm="ARK45" interval="3" tolerance="1e-3">
+			<samples>128</samples>
+			<operators>
+				<integration_vectors>density_matrix</integration_vectors>
+				<![CDATA[
+				E1c = conj(E1);
+				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
+				r11=1-r22-r33-r44; // sum of all level populations is 1
+				dr12_dt = -(g12+i*(d1-d2))*r12 - i/2.*(E2*r13+Ec*r14-E1c*r32);
+				dr13_dt = -(g13+i*d1)*r13 - i/2.*(E1c*(r11-r33)+E2c*r12);
+				dr14_dt = -(g14+i*(d1-d2+dc))*r14 - i/2.*(Ecc*r12-E1c*r34);
+				dr22_dt = -W21*r22 + W42*r44 + W32*r33 - i/2.*(Ec*r24-Ecc*r42+E2*r23-E2c*r32);
+				dr23_dt = -(g23+i*d2)*r23 - i/2.*(E2c*(r22-r33)+E1c*r21-Ecc*r43);
+				dr24_dt = -(g24+i*dc)*r24 - i/2.*(Ecc*(r22-r44)-E2c*r34);
+				dr33_dt = -(W31+W32)*r33 + W43*r44 - i/2.*(E1c*r31-E1*r13+E2c*r32-E2*r23);
+				dr34_dt = -(g34+i*(-d2+dc))*r34 - i/2.*(Ecc*r32-E1*r14-E2*r24);
+				dr44_dt = -(W41+W42+W43)*r44 - i/2.*(Ecc*r42-Ec*r24);
+				]]>
+			</operators>
+		</integrate>
+	</sequence>
+
+
+
+
+	<!-- The output to generate -->
+	<output format="binary" filename="Nlevels_no_dopler_no_z.xsil">
+		<group>
+			<sampling basis="d1" initial_sample="yes">
+				<dependencies>density_matrix</dependencies>
+				<moments>>r13_rlOut r13_imOut >r23_rlOut r23_imOut r24_rlOut r24_imOut</moments>
+				<![CDATA[
+				r13_rlOut = r13.Re();
+				r13_imOut = r13.Im();
+				r23_rlOut = r23.Re();
+				r23_imOut = r23.Im();
+				r24_rlOut = r24.Re();
+				r24_imOut = r24.Im();
+				]]>
+			</sampling>
+		</group>
+	</output>
+
+</simulation>
+	
+<!-- 
+vim: ts=2 sw=2 foldmethod=indent: 
+-->
diff --git a/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Nlevels_no_dopler_no_z_from_mathematica.xmds b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Nlevels_no_dopler_no_z_from_mathematica.xmds
new file mode 100644
index 0000000..3f306c5
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/Nlevels_no_dopler_no_z_from_mathematica.xmds
@@ -0,0 +1,178 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+	<name>Nlevels_no_dopler_no_z_from_mathematica</name>
+
+	<author>Eugeniy Mikhailov</author>
+	<description>
+		License GPL.
+
+		Solving 4 level atom in N-field configuration,
+		no field propagation along spatial axis included
+		no Doppler broadening
+
+		For master equations look "Four-level 'N-scheme' in bare and quasi-dressed states pictures"
+		by T. Abi-Salloum, S. Meiselman, J.P. Davis and F.A. Narducci
+		Journal of Modern Optics, 56: 18, 1926 -- 1932, (2009).
+
+		Present calculation matches Fig.3 (b) from the above paper. 
+		Note that I need to double all Rabi frequencies to match the figure.
+
+		*  -------- |4>
+		*    \        
+		*     \ Ec    -------- |3>  
+		*      \       /    \
+		*       \  E2 /      \ 
+		*        \   /        \ E1
+		*       ------- |2>    \
+		*                       \
+		*                    ------- |1>
+		*
+
+		We moved to dimensionless units
+		t -> t*g    ,time
+		z -> z*g/c  , distance
+		rabi_frequency -> rabi_frequency/g 
+		eta -> eta*c/g^2      , coupling constant
+		gij -> gij/g
+		Wij -> Wij/g
+
+		where g is 1MHz rate
+	</description>
+
+	<features>
+		<benchmark />
+		<bing />
+		<fftw plan="patient" />
+		<openmp />
+		<auto_vectorise />
+		<globals>
+			<![CDATA[
+			// 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;
+
+			// branching ratios
+			const double R41=0.5, R42=0.5;
+			const double R31=0.5, R32=0.5;
+
+			// const double d1=0;         // E2 detuning in MHz
+			const double d2=0;         // E2 detuning in MHz
+			const double d3=0;         // Ec detuning in MHz
+
+			const complex E1=0.005/2;     // Rabi frequency in MHz
+			const complex E2=5.0/2;         // Rabi frequency in MHz
+			const complex Ec=6.0/2;         // Rabi frequency in MHz
+
+			complex  E1c, E2c, Ecc; // Complex conjugated Rabi frequencies
+
+			complex  r21, r31, r41, r32, r42, r43, r44;
+
+			]]>
+		</globals>
+	</features>
+
+	<geometry>
+		<propagation_dimension> t </propagation_dimension>
+		<transverse_dimensions>
+			<dimension name="d1" lattice="1280"  domain="(-60, 60)" />
+		</transverse_dimensions>
+	</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</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
+				 so r11 is missing
+		-->
+		<initialisation>
+			<![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 increment 
+			// seems to be huge and adaptive solver makes smaller and smaller steps.
+			// 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
+			r11 = 1; r22 = 0; r33 = 0;
+			r12 = 0; r13 = 0; r14 = 0;
+			r23 = 0; r24 = 0;
+			r34 = 0;
+			]]>
+		</initialisation>
+	</vector>
+
+	<sequence>
+		<integrate algorithm="ARK89" interval="3" tolerance="1e-3">
+			<samples>128</samples>
+			<operators>
+				<integration_vectors>density_matrix</integration_vectors>
+				<![CDATA[
+				E1c = conj(E1);
+				E2c = conj(E2);
+				Ecc = conj(Ec);
+
+				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>
+	</sequence>
+
+
+
+
+	<!-- The output to generate -->
+	<output format="binary" filename="Nlevels_no_dopler_no_z_from_mathematica.xsil">
+		<group>
+			<sampling basis="d1" initial_sample="yes">
+				<dependencies>density_matrix</dependencies>
+				<moments>>r13_rlOut r13_imOut >r23_rlOut r23_imOut r24_rlOut r24_imOut</moments>
+				<![CDATA[
+				r13_rlOut = r13.Re();
+				r13_imOut = r13.Im();
+				r23_rlOut = r23.Re();
+				r23_imOut = r23.Im();
+				r24_rlOut = r24.Re();
+				r24_imOut = r24.Im();
+				]]>
+			</sampling>
+		</group>
+	</output>
+
+</simulation>
+	
+<!-- 
+vim: ts=2 sw=2 foldmethod=indent: 
+-->
diff --git a/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/pp.m b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/pp.m
new file mode 100644
index 0000000..ed62f63
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/pp.m
@@ -0,0 +1,23 @@
+Nlevels_no_dopler_no_z;
+
+figure(1)
+imagesc(t_1, d1_1, r13_imOut_1)
+colorbar;
+xlabel('t (uS)')
+ylabel('d_1 (MHz)')
+zlabel('r13_imOut_1')
+title('Imag(rho_{13}) vs time and d_2 detuning')
+
+
+figure(2)
+plot(d1_1, r13_imOut_1(:,end),'-')
+xlabel('d_1 (MHz)')
+ylim ([-6e-4,0])
+title('Imag(rho_{13}) vs  d_1 detuning')
+print('imag_r13.pdf');
+
+figure(3)
+plot(d1_1,r13_rlOut_1(:,end),'-')
+xlabel('d_1 (MHz)')
+title('Real(rho_{13}) vs  d_1 detuning')
+print('real_r13.pdf');
diff --git a/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/pp_from_mathematica.m b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/pp_from_mathematica.m
new file mode 100644
index 0000000..1796093
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_no_z.Frank_fig3b/pp_from_mathematica.m
@@ -0,0 +1,23 @@
+Nlevels_no_dopler_no_z_from_mathematica;
+
+figure(1)
+imagesc(t_1, d1_1, r13_imOut_1)
+colorbar;
+xlabel('t (uS)')
+ylabel('d_1 (MHz)')
+zlabel('r13_imOut_1')
+title('Imag(rho_{13}) vs time and d_2 detuning')
+
+
+figure(2)
+plot(d1_1, r13_imOut_1(:,end),'-')
+xlabel('d_1 (MHz)')
+ylim ([-6e-4,0])
+title('Imag(rho_{13}) vs  d_1 detuning')
+print('imag_r13.pdf');
+
+figure(3)
+plot(d1_1,r13_rlOut_1(:,end),'-')
+xlabel('d_1 (MHz)')
+title('Real(rho_{13}) vs  d_1 detuning')
+print('real_r13.pdf');