summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Makefile54
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Nlevels_no_dopler_with_z_4wm.xmds321
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/Makefile54
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/map2dat.m33
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/plot_fields_propagation_I2.gp15
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/pp_I2.m71
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/map2dat.m33
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp.m190
-rw-r--r--xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp_I2.m71
9 files changed, 842 insertions, 0 deletions
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Makefile
new file mode 100644
index 0000000..578f2b0
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Makefile
@@ -0,0 +1,54 @@
+### -*- 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
+ rm -f $(png_targets)
+
+eps_targets = $(wildcard *.eps)
+pdf_targets = $(eps_targets:%.eps=%.pdf)
+png_targets = $(pdf_targets:%.pdf=%.png)
+
+png: pdf $(png_targets)
+
+$(png_targets): %.png : %.pdf
+ convert -density 300 $< $@
+
+pdf: $(pdf_targets)
+
+$(pdf_targets): %.pdf : %.eps
+ cat $< | ps2eps -B > __tt.eps
+ epspdf __tt.eps $@
+ rm -f __tt.eps
+ #ps2eps -B $< | epspdf $< $@
+.PRECIOUS: %.run %.xsil %.m
+.PHONY: all clean
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Nlevels_no_dopler_with_z_4wm.xmds b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Nlevels_no_dopler_with_z_4wm.xmds
new file mode 100644
index 0000000..e0e83b3
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/Nlevels_no_dopler_with_z_4wm.xmds
@@ -0,0 +1,321 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+ <name>Nlevels_no_dopler_with_z_4wm</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="(-0.5e-6, 0.5e-6)" />
+ </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>
+ <!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
+ <!-- Comment out no light field initial conditions
+ <![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;
+ ]]>
+ -->
+ <!-- Below initialization assumes strong E1 and E3 which were shining for long time before
+ we even start to look at the problem -->
+ <!-- Precalculated by Nate via Mathematica steady state solver -->
+ <dependencies>E_field</dependencies>
+ <![CDATA[
+ E1c = conj(E1);
+ E2c = conj(E2);
+ E3c = conj(E3);
+ E4c = conj(E4);
+
+ // IMPORTANT: assumes no detunings
+ r11 = (gt*(4*mod2(E1) + (G3 + gt)*(G3 + 2*gt))*((G4 +
+ 2*gt)*(4*mod2(E3) + gt*(G4 + gt)) + 4*mod2(E3)*G4*(R41 -
+ R42)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
+ 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
+ gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
+ 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
+
+ r13 = -((E1c*gt*(G3 + gt)*i*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt))
+ + 4*mod2(E3)*G4*(R41 -
+ R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
+ 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
+ gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
+ 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
+
+ r22 = (gt*(4*mod2(E3) + (G4 + gt)*(G4 + 2*gt))*((G3 +
+ 2*gt)*(4*mod2(E1) + gt*(G3 + gt)) + 4*mod2(E1)*G3*(-R31 +
+ R32)))/(2.*(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
+ 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
+ gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
+ 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
+
+ r24 = -((E3c*gt*(G4 + gt)*i*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt))
+ + 4*mod2(E1)*G3*(-R31 +
+ R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
+ 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
+ gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
+ 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42)));
+
+ r33 = (2*mod2(E1)*gt*((G4 + 2*gt)*(4*mod2(E3) + gt*(G4 + gt)) +
+ 4*mod2(E3)*G4*(R41 -
+ R42)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
+ 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
+ gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
+ 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42));
+
+ r44 = (2*mod2(E3)*gt*((G3 + 2*gt)*(4*mod2(E1) + gt*(G3 + gt)) +
+ 4*mod2(E1)*G3*(-R31 +
+ R32)))/(-16*mod2(E1)*mod2(E3)*G3*G4*R32*R41 + (G3 +
+ 2*gt)*(4*mod2(E1) + gt*(G3 + gt))*((G4 + 2*gt)*(4*mod2(E3) +
+ gt*(G4 + gt)) - 4*mod2(E3)*G4*R42) + 4*mod2(E1)*G3*R31*(-((G4 +
+ 2*gt)*(4*mod2(E3) + gt*(G4 + gt))) + 4*mod2(E3)*G4*R42));
+
+ ]]>
+ </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">
+ <!--This set boundary condition at all 'z' and left border of 't' (i.e. min(t))-->
+ <!--
+ <![CDATA[
+ r11 = 0; r22 = 1; r33 = 0; r44 = 0;
+ r12 = 0; r13 = 0; r14 = 0;
+ r23 = 0; r24 = 0;
+ r34 = 0;
+ printf("z= %g, t= %g\n", z, t);
+ ]]>
+ -->
+ </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_4wm.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_for_irina_pqe/fast_light/Makefile b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/Makefile
new file mode 100644
index 0000000..76d3f5c
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/Makefile
@@ -0,0 +1,54 @@
+### -*- make -*-
+### This makefile can be used to build and run the XMDS examples
+
+
+XSIL_FILES = Nlevels_no_dopler_with_z_4wm.xsil
+M_FILES = $(patsubst %.xsil,%.m,$(XSIL_FILES))
+GNUPLOT_FILES = $(wildcard *.gp)
+
+XSIL2GRAPHICS = xsil2graphics
+
+# fast light
+PARAMS = --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=3.8e7 --E4o=6.3e4
+# slow light EIT
+#PARAMS = --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=0 --E4o=0
+
+all: $(XSIL_FILES) Nlevels_no_dopler_with_z_4wm.xsil $(M_FILES) plot png
+
+Nlevels_no_dopler_with_z_4wm.xsil: ../Nlevels_no_dopler_with_z_4wm.run
+ $< $(PARAMS) | grep "Time elapsed for simulation is:" > exact_analysis_execution_time.txt
+
+%.m: %.xsil
+ $(XSIL2GRAPHICS) $<
+
+plot: $(M_FILES) $(GNUPLOT_FILES)
+ octave pp_I2.m
+ gnuplot plot_fields_propagation_I2.gp
+
+clean:
+ rm -f $(CC_FILES) $(RUN_FILES) $(M_FILES) $(XSIL_FILES) *.wisdom.fftw3 *.dat octave-core *.wisdom *.pdf
+ rm -f $(pdf_targets)
+ rm -f $(eps_targets)
+
+real_clean: clean
+ rm -f $(png_targets)
+
+eps_targets = $(wildcard *.eps)
+pdf_targets = $(eps_targets:%.eps=%.pdf)
+png_targets = $(pdf_targets:%.pdf=%.png)
+
+pdf: $(pdf_targets)
+
+$(pdf_targets): %.pdf : %.eps
+ cat $< | ps2eps -B > __tt.eps
+ epspdf __tt.eps $@
+ rm -f __tt.eps
+ #ps2eps -B $< | epspdf $< $@
+
+png: pdf $(png_targets)
+
+$(png_targets): %.png : %.pdf
+ convert -density 300 $< $@
+
+.PRECIOUS: %.run %.xsil %.m
+.PHONY: all clean
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/map2dat.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/map2dat.m
new file mode 100644
index 0000000..969b6dc
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/map2dat.m
@@ -0,0 +1,33 @@
+function map2dat(outfile, x,y,z, xskip, yskip)
+% saves 3D data in suitable way to be drawn by gnuplot
+% x,y - vectors of x,y values
+% z map of z values as used by Octave/Matlab
+% xskip, yskip - skip paprameters
+% only every, xskip, yskip point will be written
+
+
+
+Nx=length(x);
+Ny=length(y);
+Nxs=Nx/xskip;
+Nys=Ny/yskip;
+points=zeros(1,3*Nxs*Nys);
+%points=[];
+tic;
+for i=1:Nxs
+ for k=1:Nys
+ %points=[points x(i*xskip) y(k*yskip) z(k*yskip,i*xskip)];
+ points((i-1)*(Nys-1)*3+3*(k-1)+1) = x(i*xskip);
+ points((i-1)*(Nys-1)*3+3*(k-1)+2) = y(k*yskip);
+ points((i-1)*(Nys-1)*3+3*(k-1)+3) = z(k*yskip,i*xskip);
+ end
+end
+disp('=== points formation complete ===');
+toc
+tic;
+%points
+fd = fopen(outfile, "wt");
+fprintf (fd, "%g %g %g\n", points);
+fclose(fd);
+disp('=== points saving complete ===');
+toc;
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/plot_fields_propagation_I2.gp b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/plot_fields_propagation_I2.gp
new file mode 100644
index 0000000..8591d87
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/plot_fields_propagation_I2.gp
@@ -0,0 +1,15 @@
+set terminal postscript portrait enhanced color solid size 5,3.5
+set output 'fields_propagation_I2.eps'
+set dgrid3d 100,100 qnorm 4
+set pm3d map
+#set contour
+set hidden3d
+set palette rgb 10,13,31 negative
+
+
+set xlabel "z (cm)"
+set ylabel "t ({/Symbol m}S)"
+set zlabel "I_2 (1/S)"
+set nokey
+#set view map
+splot [0:1.5][-0.2:0.2] 'I2.dat'
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/pp_I2.m
new file mode 100644
index 0000000..cb34777
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/fast_light/pp_I2.m
@@ -0,0 +1,71 @@
+Nlevels_no_dopler_with_z_4wm
+
+%% field propagation
+z_1=z_1*100; % z in cm
+t_1=t_1*1e6; % time now measured in uS
+figure(1)
+%set(gca,'fontsize',20);
+imagesc(z_1, t_1, I2_out_1); colorbar
+tmin=-0.4;
+tmax= 0.4;
+ylim([tmin,tmax],'manual');
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_2')
+title('I_2')
+
+xskip=1;
+yskip=10;
+map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
+
+
+
+print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I2.eps')
+
+
+
+%% fields before and after the cell
+figure(2)
+%set(gca,'fontsize',30);
+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')
+legend('location', 'southwest');
+
+[b, a]=butter(3, 0.05);
+I2_out_after=I2_out_1(:,end);
+I2_out_after_filtered=filtfilt(b,a,I2_out_after);
+settling_time=0.01; %uS
+t_good_indx=t_1> min(t_1 + settling_time);
+[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
+delay_time=t_1(max_pos_after)-t_1(max_pos_before);
+printf('Second field delay time = %f uS\n',delay_time);
+
+%set(gca,'fontsize',40);
+%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
+print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell_I2.eps')
+
+figure(4)
+I2_max_in=max(I2_out_1(t_good_indx,1));
+I2_max_out=max(I2_out_1(t_good_indx,end));
+I2_in_norm=(I2_out_1(:,1))/I2_max_in;
+I2_out_norm=(I2_out_1(:,end))/I2_max_out;
+tmin=-.05;
+tmax=.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 ...
+ )
+legend('location', 'southeast');
+xlim([tmin,tmax],'manual');
+xlabel('t (uS)')
+ylabel('I_2')
+title('I_2 before and after cell normalized')
+%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
+print('-color','-depsc2', '-tight','-S200,120', 'probe_before_after_cell_I2_normalized.eps')
+
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/map2dat.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/map2dat.m
new file mode 100644
index 0000000..969b6dc
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/map2dat.m
@@ -0,0 +1,33 @@
+function map2dat(outfile, x,y,z, xskip, yskip)
+% saves 3D data in suitable way to be drawn by gnuplot
+% x,y - vectors of x,y values
+% z map of z values as used by Octave/Matlab
+% xskip, yskip - skip paprameters
+% only every, xskip, yskip point will be written
+
+
+
+Nx=length(x);
+Ny=length(y);
+Nxs=Nx/xskip;
+Nys=Ny/yskip;
+points=zeros(1,3*Nxs*Nys);
+%points=[];
+tic;
+for i=1:Nxs
+ for k=1:Nys
+ %points=[points x(i*xskip) y(k*yskip) z(k*yskip,i*xskip)];
+ points((i-1)*(Nys-1)*3+3*(k-1)+1) = x(i*xskip);
+ points((i-1)*(Nys-1)*3+3*(k-1)+2) = y(k*yskip);
+ points((i-1)*(Nys-1)*3+3*(k-1)+3) = z(k*yskip,i*xskip);
+ end
+end
+disp('=== points formation complete ===');
+toc
+tic;
+%points
+fd = fopen(outfile, "wt");
+fprintf (fd, "%g %g %g\n", points);
+fclose(fd);
+disp('=== points saving complete ===');
+toc;
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp.m
new file mode 100644
index 0000000..b88fa6e
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp.m
@@ -0,0 +1,190 @@
+Nlevels_no_dopler_with_z_4wm
+
+%% field propagation
+z_1=z_1*100; % z in cm
+t_1=t_1*1e6; % time now measured in uS
+figure(1)
+subplot(2,2,1); imagesc(z_1, t_1, I1_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_1')
+title('I_1')
+subplot(2,2,2); imagesc(z_1, t_1, I2_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_2')
+title('I_2')
+subplot(2,2,3); imagesc(z_1, t_1, I3_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_3')
+title('I_3')
+subplot(2,2,4); imagesc(z_1, t_1, I4_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_4')
+title('I_4')
+
+
+print('-color','fields_propagation.eps')
+
+
+
+%% fields before and after the cell
+figure(2)
+subplot(2,2,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(2,2,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(2,2,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);
+settling_time=0.8; %uS
+t_good_indx=t_1> min(t_1 + settling_time);
+[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
+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')
+
+subplot(2,2,4);
+plot( ...
+ t_1,I4_out_1(:,1),'-;before;', "linewidth", 4, ...
+ t_1,I4_out_1(:,end), '-;after;', "linewidth", 4 ...
+ )
+xlabel('t (uS)')
+ylabel('I_3 (1/s)^2')
+title('I_3 before and after cell')
+
+figure(4)
+I2_max_in=max(I2_out_1(t_good_indx,1));
+I2_max_out=max(I2_out_1(t_good_indx,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})')
+
+
diff --git a/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp_I2.m b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp_I2.m
new file mode 100644
index 0000000..45d4913
--- /dev/null
+++ b/xmds2/Nlevels_no_dopler_with_z_4wm_for_irina_pqe/pp_I2.m
@@ -0,0 +1,71 @@
+Nlevels_no_dopler_with_z_4wm
+
+%% field propagation
+z_1=z_1*100; % z in cm
+t_1=t_1*1e6; % time now measured in uS
+figure(1)
+%set(gca,'fontsize',20);
+imagesc(z_1, t_1, I2_out_1); colorbar
+tmin=-0.4;
+tmax= 0.4;
+ylim([tmin,tmax],'manual');
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_2')
+title('I_2')
+
+xskip=1;
+yskip=10;
+%map2dat('I2.dat',z_1,t_1, I2_out_1, xskip, yskip);
+
+
+
+print('-color','-depsc2', '-tight', '-S200,120', 'fields_propagation_I2.eps')
+
+
+
+%% fields before and after the cell
+figure(2)
+%set(gca,'fontsize',30);
+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')
+legend('location', 'southwest');
+
+[b, a]=butter(3, 0.05);
+I2_out_after=I2_out_1(:,end);
+I2_out_after_filtered=filtfilt(b,a,I2_out_after);
+settling_time=0.8; %uS
+t_good_indx=t_1> min(t_1 + settling_time);
+[m,max_pos_before]=max(I2_out_1(t_good_indx,1) ); [m,max_pos_after]=max(I2_out_after_filtered(t_good_indx));
+delay_time=t_1(max_pos_after)-t_1(max_pos_before);
+printf('Second field delay time = %f uS\n',delay_time);
+
+%set(gca,'fontsize',40);
+%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
+print('-color','-depsc2', '-tight','-S200,120', 'fields_before_after_cell_I2.eps')
+
+figure(4)
+I2_max_in=max(I2_out_1(t_good_indx,1));
+I2_max_out=max(I2_out_1(t_good_indx,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 ...
+ )
+legend('location', 'southeast');
+xlim([tmin,tmax],'manual');
+xlabel('t (uS)')
+ylabel('I_2')
+title('I_2 before and after cell normalized')
+%set (gcf,'paperposition',[0.5 0 2.5,1.5]); % IMPORTANT to shrink eps size for readable fonts
+print('-color','-depsc2', '-tight','-S200,120', 'probe_before_after_cell_I2_normalized.eps')
+