four-wave mixing added

3 files changed, 469 insertions, 0 deletions

diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm/Makefile
new file mode 100644
index 0000000..a89800a
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/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_4wm/Nlevels_no_dopler_with_z.xmds b/xmds2/Nlevels_no_dopler_with_z_4wm/Nlevels_no_dopler_with_z.xmds
new file mode 100644
index 0000000..b17ef61
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/Nlevels_no_dopler_with_z.xmds
@@ -0,0 +1,259 @@
+<?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
+
+		We assume four-wave mixing condition when w3-w4=w2-w1 i.e. fields E3 and E4 drive the same
+		resonance as fields E2 and E1.
+
+
+		*  --------------- |4>
+		*    \        \    
+		*     \ E3     \    -------- |3>  
+		*      \     E4 \   /    \
+		*       \        \ / E2   \ 
+		*        \        /        \ E1
+		*   |2> --------------      \
+		*                   \        \
+		*                    \        \
+		*                   ------------- |1>
+		*
+
+
+		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=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)]
+
+				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
+				const double gt=0.01 *(2*M_PI*1e6);
+				// Natural linewidth  of j's level in [1/s]
+				const double G3=3.0 *(2*M_PI*1e6);
+				const double G4=3.0 *(2*M_PI*1e6);
+
+				// 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;
+
+
+				complex  E1c, E2c, E3c, E4c; // Complex conjugated Rabi frequencies
+
+				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
+			]]>
+		</globals>
+		<benchmark />
+    <arguments>
+			<!-- Rabi frequency divided by 2 in [1/s] -->
+      <argument name="E1o" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
+      <argument name="E2o" type="real" default_value="0.05*(2*M_PI*1e6)" />
+      <argument name="E3o" type="real" default_value="2*3.0*(2*M_PI*1e6)" />
+      <argument name="E4o" type="real" default_value=".01*(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" />
+    </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="(-14.0e-7, 4.0e-7)" />
+		</transverse_dimensions>
+	</geometry>
+
+	<!-- Rabi frequency --> 
+	<vector name="E_field" type="complex" initial_space="t">
+		<components>E1 E2 E3 E4</components>
+		<initialisation>
+			<![CDATA[
+			// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
+			// as well as time
+			E1=E1o;
+			E2=E2o*exp(-pow( ((t-0.0)/.1e-6),2) );
+			E3=E3o;
+			E4=E4o;
+			]]>
+		</initialisation>
+	</vector>
+
+	<vector name="density_matrix" type="complex" initial_space="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 r44</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; r44 = 0;
+			// r12 = 0; r13 = 0; r14 = 0;
+			// r23 = 0; r24 = 0;
+			// r34 = 0;
+			]]>
+		</initialisation>
+	</vector>
+
+	<sequence>
+		<!--For this set of conditions ARK45 is faster than ARK89-->
+		<integrate algorithm="ARK45" tolerance="1e-5" interval="1.5e-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-->
+
+		<!--<integrate algorithm="SIC" interval="4e-2" steps="200">-->
+			<samples>200 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 = 0; r22 = 1; r33 = 0; r44 = 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);
+						E4c = conj(E4);
+
+						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 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 + delta2)*r12 - E2*r13 - E3*r14 + E1c*r32 + E4c*r42);
+						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(E1c*r11) - E2c*r12 - delta1*r13 + E1c*r33 + E4c*r43);
+						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-(E4c*r11) - E3c*r12 - (delta1 - delta2 + delta3)*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*r23 + E2c*r33 + E3c*r43);
+						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(E4c*r21) - E3c*r22 - delta3*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 - delta3)*r34);
+						dr44_dt = i*(E4*r14 + E3*r24 - E4c*r41 - E3c*r42) - (G4 + 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[
+					dE1_dz = i*eta*conj(r13) -Lt[E1] ;
+					dE2_dz = i*eta*conj(r23) -Lt[E2] ;
+					dE3_dz = i*eta*conj(r24) -Lt[E3] ;
+					dE4_dz = i*eta*conj(r14) -Lt[E4] ;
+          ]]>
+			</operators>
+		</integrate>
+	</sequence>
+
+
+
+
+	<!-- The output to generate -->
+	<output format="binary" filename="Nlevels_no_dopler_with_z.xsil">
+		<group>
+      <sampling basis="t(1000)" initial_sample="yes">
+				<dependencies>E_field</dependencies>
+				<moments>I1_out I2_out I3_out I4_out</moments>
+				<![CDATA[
+				I1_out = mod2(E1);
+				I2_out = mod2(E2);
+				I3_out = mod2(E3);
+				I4_out = mod2(E4);
+				]]>
+			</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: 
+-->
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm/pp.m b/xmds2/Nlevels_no_dopler_with_z_4wm/pp.m
new file mode 100644
index 0000000..ab536be
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm/pp.m
@@ -0,0 +1,173 @@
+Nlevels_no_dopler_with_z
+
+%% field propagation
+z_1=z_1*100; % z in cm
+t_1=t_1*1e6; % time now measured in uS
+figure(1)
+subplot(1,3,1); imagesc(z_1, t_1, I1_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_1')
+title('I_1')
+subplot(1,3,2); imagesc(z_1, t_1, I2_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_2')
+title('I_2')
+subplot(1,3,3); imagesc(z_1, t_1, I3_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_3')
+title('I_3')
+
+print('-color','fields_propagation.eps')
+
+
+
+%% fields before and after the cell
+figure(2)
+subplot(1,3,1);
+plot( ...
+	t_1,I1_out_1(:,1),'-;before;', "linewidth", 4,  ...
+	t_1,I1_out_1(:,end), '-;after;', "linewidth", 4 ...
+	)
+xlabel('t (uS)')
+ylabel('I_1 (1/s)^2')
+title('I_1 before and after cell')
+subplot(1,3,2);
+plot( ...
+	t_1,I2_out_1(:,1),'-;before;', "linewidth", 4,  ...
+	t_1,I2_out_1(:,end), '-;after;', "linewidth", 4 ...
+	)
+xlabel('t (uS)')
+ylabel('I_2 (1/s)^2')
+title('I_2 before and after cell')
+subplot(1,3,3);
+plot( ...
+	t_1,I3_out_1(:,1),'-;before;', "linewidth", 4,  ...
+	t_1,I3_out_1(:,end), '-;after;', "linewidth", 4 ...
+	)
+xlabel('t (uS)')
+ylabel('I_3 (1/s)^2')
+title('I_3 before and after cell')
+
+[b, a]=butter(3, 0.05);
+I2_out_after=I2_out_1(:,end);
+I2_out_after_filtered=filtfilt(b,a,I2_out_after);
+[m,max_pos_before]=max(I2_out_1(:,1) ); [m,max_pos_after]=max(I2_out_after_filtered);
+delay_time=t_1(max_pos_after)-t_1(max_pos_before);
+printf('Second field delay time  = %f uS\n',delay_time);
+
+print('-color','fields_before_after_cell.eps')
+
+figure(4)
+I2_max_in=max(I2_out_1(:,1));
+I2_max_out=max(I2_out_1(:,end));
+I2_in_norm=(I2_out_1(:,1))/I2_max_in;
+I2_out_norm=(I2_out_1(:,end))/I2_max_out;
+tmin=-0.05;
+tmax=0.05;
+indx=(t_1>=tmin & t_1<=tmax);  % soom in in time to this region
+plot( ...
+	t_1(indx),I2_in_norm(indx),'-;before;', "linewidth", 4,  ...
+	t_1(indx),I2_out_norm(indx), '-;after;', "linewidth", 4 ...
+	)
+xlim([tmin,tmax],'manual');	
+xlabel('t (uS)')
+ylabel('I_2')
+title('I_2 before and after cell normalized')
+print('-color','probe_before_after_cell.eps')
+
+%return;
+
+%% all density matrix elements in one plot
+% diagonal populations, 
+% upper triangle real part of coherences, 
+% lower diagonal imaginary part of coherences
+z_2=z_2*100; % z in cm
+t_2=t_2*1e6; % time now measured in uS
+figure(3)
+subplot(4,4,1); imagesc (z_2, t_2, r11_out_2);  caxis([0,1]); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('rho_{11}')
+title('rho_{11}')
+subplot(4,4,6); imagesc (z_2, t_2, r22_out_2);  caxis([0,1]); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('rho_{22}')
+title('rho_{22}')
+subplot(4,4,11); imagesc (z_2, t_2, r33_out_2);  caxis([0,1]); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('rho_{33}')
+title('rho_{33}')
+subplot(4,4,16); imagesc (z_2, t_2, r44_out_2); caxis([0,1]);  colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('rho_{44}')
+title('rho_{44}')
+% real parts of coherences
+subplot(4,4,2); imagesc(z_2, t_2, r12_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{12})')
+title('Real(rho_{12})')
+subplot(4,4,3); imagesc(z_2, t_2, r13_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{13})')
+title('Real(rho_{13})')
+subplot(4,4,4); imagesc(z_2, t_2, r14_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{14})')
+title('Real(rho_{14})')
+subplot(4,4,7); imagesc(z_2, t_2, r23_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{23})')
+title('Real(rho_{23})')
+subplot(4,4,8); imagesc(z_2, t_2, r24_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{24})')
+title('Real(rho_{24})')
+subplot(4,4,12); imagesc(z_2, t_2, r34_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{34})')
+title('Real(rho_{34})')
+% imaginary parts of coherences
+subplot(4,4,5); imagesc(z_2, t_2, r12_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{12})')
+title('Imag(rho_{12})')
+subplot(4,4,9); imagesc(z_2, t_2, r13_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{13})')
+title('Imag(rho_{13})')
+subplot(4,4,10); imagesc(z_2, t_2, r23_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{23})')
+title('Imag(rho_{23})')
+subplot(4,4,13); imagesc(z_2, t_2, r14_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{14})')
+title('Imag(rho_{14})')
+subplot(4,4,14); imagesc(z_2, t_2, r24_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{24})')
+title('Imag(rho_{24})')
+subplot(4,4,15); imagesc(z_2, t_2, r34_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{34})')
+title('Imag(rho_{34})')
+
+