summaryrefslogtreecommitdiff
path: root/xmds2
diff options
context:
space:
mode:
Diffstat (limited to 'xmds2')
-rw-r--r--xmds2/Nlevels_no_dopler_with_z/Makefile37
-rw-r--r--xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds228
-rw-r--r--xmds2/Nlevels_no_dopler_with_z/pp.m8
3 files changed, 273 insertions, 0 deletions
diff --git a/xmds2/Nlevels_no_dopler_with_z/Makefile b/xmds2/Nlevels_no_dopler_with_z/Makefile
new file mode 100644
index 0000000..a89800a
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z/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_with_z/Nlevels_no_dopler_with_z.xmds b/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds
new file mode 100644
index 0000000..fc4da89
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z/Nlevels_no_dopler_with_z.xmds
@@ -0,0 +1,228 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+ <name>Nlevels_no_dopler_with_z</name>
+
+ <author>Eugeniy Mikhailov</author>
+ <description>
+ License GPL.
+
+ Solving 4 level atom in N-field configuration,
+ with field propagation along spatial axis Z
+ 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>
+ * \
+ * \ E3 -------- |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
+
+ We are solving
+ dE/dz+(1/c)*dE/dt=i*eta*rho_ij, where j llevel is higher then i
+ in xmds terms it looks like
+ dE_dz = i*eta*rhoij - L[E], here we moved t dependence to Fourier space
+ </description>
+
+ <features>
+ <benchmark />
+ <arguments>
+ <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" />
+ <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" />
+ </arguments>
+ <bing />
+ <fftw plan="patient" />
+ <openmp />
+ <auto_vectorise />
+ <globals>
+ <![CDATA[
+ const complex UpperDecayRate=10e6; // measured in s^-1
+ const double c=3e8;
+ const double lambda=794.7e-9; //wavelength in m
+ const double N=1e9; //number of particles per cubic cm
+ const complex eta = 3*lambda*lambda*(N/1e-6)*UpperDecayRate/8.0/3.14*(c/UpperDecayRate/UpperDecayRate); //dimensionless coupling constant
+
+ // 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 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;
+ ]]>
+ </globals>
+ </features>
+
+ <geometry>
+ <propagation_dimension> z </propagation_dimension>
+ <transverse_dimensions>
+ <dimension name="t" lattice="128" domain="(-60, 60)" />
+ </transverse_dimensions>
+ </geometry>
+
+ <vector name="E_field" type="complex" dimensions="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;
+ E2=E2o;
+ E3=E3o;
+ ]]>
+ </initialisation>
+ </vector>
+
+ <vector name="density_matrix" type="complex" dimensions="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>-->
+ <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="SIC" interval="4" steps="200">
+ <samples>200</samples>
+ <operators>
+ <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
+ <integration_vectors>density_matrix</integration_vectors>
+ <dependencies>E_field</dependencies>
+ <boundary_condition kind="left">
+ <![CDATA[
+ r11 = 1; r22 = 0; r33 = 0;
+ r12 = 0; r13 = 0; r14 = 0;
+ r23 = 0; r24 = 0;
+ r34 = 0;
+ ]]>
+ </boundary_condition>
+ <![CDATA[
+ E1c = conj(E1);
+ E2c = conj(E2);
+ E3c = conj(E3);
+
+ 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*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1*r32);
+ dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1*r11) - E2*r12 - delta1*r13 + E1*r33);
+ dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E3*r12) - (delta1 - delta2 + delta3)*r14 + E1*r34);
+ //dr21_dt = -(gt*r21) + i*((delta1 - delta2)*r21 - E1*r23 + E2*r31 + E3*r41);
+ dr22_dt = gt/2. - gt*r22 + i*(-(E2*r23) - E3*r24 + E2*r32 + E3*r42) + G3*r33*R32 + G4*r44*R42;
+ dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(E1*r21) - E2*r22 - delta2*r23 + E2*r33 + E3*r43);
+ dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E3*r22) - delta3*r24 + E2*r34 + E3*r44);
+ //dr31_dt = -((G3 + 2*gt)*r31)/2. + i*(E1*r11 + E2*r21 + delta1*r31 - E1*r33);
+ //dr32_dt = -((G3 + 2*gt)*r32)/2. + i*(E1*r12 + E2*r22 + delta2*r32 - E2*r33 - E3*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 - E3*r32 + (delta2 - delta3)*r34);
+ //dr41_dt = -((G4 + 2*gt)*r41)/2. + i*(E3*r21 + (delta1 - delta2 + delta3)*r41 - E1*r43);
+ //dr42_dt = -((G4 + 2*gt)*r42)/2. + i*(E3*r22 + delta3*r42 - E2*r43 - E3*r44);
+ //dr43_dt = -((G3 + G4 + 2*gt)*r43)/2. + i*(E3*r23 - E1*r41 - E2*r42 + (-delta2 + delta3)*r43);
+ //dr44_dt = E3*i*(r24 - r42) - (G4 + gt)*r44;
+ ]]>
+ </operator>
+ <operator kind="ex" constant="yes">
+ <operator_names>Lt</operator_names>
+ <![CDATA[
+ Lt = i*kt;
+ ]]>
+ </operator>
+ <integration_vectors>E_field</integration_vectors>
+ <![CDATA[
+ dE1_dz = i*eta*r13 - Lt[E1];
+ dE2_dz = i*eta*r23 - Lt[E2];
+ dE3_dz = i*eta*r24 - Lt[E3];
+ ]]>
+ </operators>
+ </integrate>
+ </sequence>
+
+
+
+
+ <!-- The output to generate -->
+ <output format="binary" filename="Nlevels_no_dopler_with_z.xsil">
+ <group>
+ <sampling basis="t" initial_sample="yes">
+ <dependencies>E_field</dependencies>
+ <moments>>I1_out I2_out I3_out</moments>
+ <![CDATA[
+ I1_out = pow(abs(E1/E1o),2);
+ I2_out = pow(abs(E2/E2o),2);
+ I3_out = pow(abs(E3/E3o),2);
+ ]]>
+ </sampling>
+ </group>
+ </output>
+
+</simulation>
+
+<!--
+vim: ts=2 sw=2 foldmethod=indent:
+-->
diff --git a/xmds2/Nlevels_no_dopler_with_z/pp.m b/xmds2/Nlevels_no_dopler_with_z/pp.m
new file mode 100644
index 0000000..30d4736
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z/pp.m
@@ -0,0 +1,8 @@
+Nlevels_no_dopler_with_z
+
+figure(1)
+imagesc(z_1, t_1, I1_out_1)
+xlabel('z')
+ylabel('t')
+zlabel('I1_{out}')
+title('I1_{out}')