Shahriar_system xmds is ready

Still need to find good set of parameters to show fast light

Programmed code includes: xmds file, pp.m to show results and some
modifications of Simon's mathematica code to better suit xmds. Total
work time 1h30m including some EIT test cases.

4 files changed, 477 insertions, 72 deletions

diff --git a/xmds2/Shahriar_system/Nlevels_no_dopler_with_z.xsil b/xmds2/Shahriar_system/Nlevels_no_dopler_with_z.xsil
new file mode 100644
index 0000000..a3b05f2
--- /dev/null
+++ b/xmds2/Shahriar_system/Nlevels_no_dopler_with_z.xsil
@@ -0,0 +1,298 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+        <name>Shahriar_system</name>
+
+        <author>Eugeniy Mikhailov, Simon Rochester</author>
+        <description>
+                License GPL.
+
+                Solving 3 level atom in double drive configuration
+                after Shahriar paper about white cavity
+                with field propagation along spatial axis Z
+                no Doppler broadening.
+
+                All fields detuned from upper level i.e. Raman configuration
+
+
+                *
+                *         .....
+                *           /   ....
+                *          / .... \
+                *         /   /    \
+                *        /   /-------- |3>  
+                *       E3  /        \
+                *      /   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=1e10*(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 *(2*M_PI*1e6);
+                                // Natural linewidth  of the upper level [1/s]
+                                const double G=2*6 *(2*M_PI*1e6);
+
+                                // incoherent pumping rate from level |1> to |3> in [1/s]
+                                const double gp=2*.6 *(2*M_PI*1e6);
+
+                                // total decay of i-th level branching ratios. Rij branching of i-th level to j-th
+                                const double R31=0.5, R32=0.5;
+
+
+                                complex  E1c, E2c, E3c; // Complex conjugated Rabi frequencies
+
+                                complex  r21, r31, r32;  // density matrix elements
+                        ]]>
+                </globals>
+                <benchmark />
+    <arguments>
+                        <!-- Rabi frequency divided by 2 in [1/s] -->
+                        <!--probe-->
+      <argument name="E1o" type="real" default_value="0.0025*(2*M_PI*1e6)" />
+                        <!--pump fields-->
+      <argument name="E2o" type="real" default_value="6.0*(2*M_PI*1e6)" />
+      <argument name="E3o" type="real" default_value="0.0*(2*M_PI*1e6)" />
+                        <!-- Fields detuning in [1/s] -->
+                        <!-- probe field detuning-->
+      <argument name="d1"  type="real" default_value="0*(2*M_PI*1e6)" />
+                        <!-- averaged detuning of pump fields i.e. mid point -->
+      <argument name="da"  type="real" default_value="0*(2*M_PI*1e6)" />
+                        <!-- detuning of pump fields with respect to each other -->
+      <argument name="delta"  type="real" default_value=".001*(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="(-2.0e-6, 4.0e-6)" />
+                </transverse_dimensions>
+        </geometry>
+
+        <!-- Rabi frequency --> 
+        <vector name="E_field" type="complex" initial_space="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*exp(-pow( ((t-0.0)/1e-6),2) );
+                        E2=E2o;
+                        E3=E3o;
+                        ]]>
+                </initialisation>
+        </vector>
+
+        <vector name="density_matrix" type="complex" initial_space="t">
+                <components>r11 r22 r33 r12 r13  r23 </components>
+                <!--
+                                 note one of the level population is redundant since
+                                 r11+r22+r33=1
+                -->
+                <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;
+                        r23 = 0;
+                        ]]>
+                </initialisation>
+        </vector>
+
+        <vector name="pump_detunings" type="complex" initial_space="t">
+                <components> d dc </components>
+                <!--dc is probably redundant since it is just complex conjugate of d-->
+                <initialisation>
+                  <![CDATA[
+                        d  = exp( i*t*delta );
+                        dc = conj(d);
+                  ]]>
+                </initialisation>
+        </vector>       
+
+
+        <sequence>
+                <!--For this set of conditions ARK45 is faster than ARK89-->
+                <integrate algorithm="ARK45" tolerance="1e-5" interval="7e-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 pump_detunings</dependencies>
+          <boundary_condition kind="left">
+            <![CDATA[
+                                                        r11 = 1; r22 = 0; r33 = 0; 
+                                                        r12 = 0; r13 = 0; 
+                                                        r23 = 0; 
+            ]]>
+          </boundary_condition>
+                                        <![CDATA[
+                                                E1c = conj(E1);
+                                                E2c = conj(E2);
+                                                E3c = conj(E3);
+
+                                                r21=conj(r12);
+                                                r31=conj(r13);
+                                                r32=conj(r23);
+
+                                                // Equations of motions according to Simon's mathematica code
+                                                dr11_dt = gt - 2*(gp + gt)*r11 - E1*i*r13 + E1c*i*r31 + G*r33;
+                                                dr12_dt = (-gp - 2*gt - d1*i + da*i)*r12 - i*(E2*d + E3*dc)*r13 + E1c*i*r32;
+                                                dr13_dt = -(E1c*i*r11) - i*(E3c*d + E2c*dc)*r12 + (-G - gp - 2*gt - d1*i)*r13 + E1c*i*r33;
+                                                dr22_dt = gt - 2*gt*r22 - i*(E2*d + E3*dc)*r23 + i*(E3c*d + E2c*dc)*r32 + G*r33;
+                                                dr23_dt = -(E1c*i*r21) - i*(E3c*d + E2c*dc)*r22 + (-G - 2*gt - da*i)*r23 + i*(E3c*d + E2c*dc)*r33;
+                                                dr33_dt = 2*gp*r11 + E1*i*r13 + i*(E2*d + E3*dc)*r23 - E1c*i*r31 - i*(E3c*d + E2c*dc)*r32 - 2*(G + gt)*r33;
+
+                                        ]]>
+        </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(r23) -Lt[E3] ;
+          ]]>
+                        </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</moments>
+                                <![CDATA[
+                                I1_out = mod2(E1);
+                                I2_out = mod2(E2);
+                                I3_out = mod2(E3);
+                                ]]>
+                        </sampling>
+                </group>
+
+                <group>
+      <sampling basis="t(100)" initial_sample="yes">
+                                <dependencies>density_matrix</dependencies>
+                                <moments>
+                                        r11_out r22_out r33_out
+                                        r12_re_out r12_im_out r13_re_out r13_im_out
+                                                              r23_re_out r23_im_out
+                                </moments>
+                                <![CDATA[
+                                // populations output 
+                                r11_out = r11.Re();
+                                r22_out = r22.Re();
+                                r33_out = r33.Re();
+                                // coherences output 
+                                r12_re_out = r12.Re();
+                                r12_im_out = r12.Im();
+                                r13_re_out = r13.Re();
+                                r13_im_out = r13.Im();
+                                r23_re_out = r23.Re();
+                                r23_im_out = r23.Im();
+                                ]]>
+                        </sampling>
+                </group>
+        </output>
+
+
+<info>
+Script compiled with XMDS2 version 2.0 "Shiny!" (HEAD)
+See http://www.xmds.org for more information.
+
+Variables that can be specified on the command line:
+  Command line argument E1o = 1.570796e+04
+  Command line argument E2o = 3.769911e+07
+  Command line argument E3o = 0.000000e+00
+  Command line argument d1 = 0.000000e+00
+  Command line argument da = 0.000000e+00
+  Command line argument delta = 6.283185e+03
+</info>
+
+<XSIL Name="moment_group_1">
+  <Param Name="n_independent">2</Param>
+  <Array Name="variables" Type="Text">
+    <Dim>5</Dim>
+    <Stream><Metalink Format="Text" Delimiter=" \n"/>
+z t I1_out I2_out I3_out 
+    </Stream>
+  </Array>
+  <Array Name="data" Type="double">
+    <Dim>201</Dim>
+    <Dim>1000</Dim>
+    <Dim>5</Dim>
+    <Stream><Metalink Format="Binary" UnsignedLong="uint32" precision="double" Type="Remote" Encoding="LittleEndian"/>
+Nlevels_no_dopler_with_z_mg0.dat
+    </Stream>
+  </Array>
+</XSIL>
+
+<XSIL Name="moment_group_2">
+  <Param Name="n_independent">2</Param>
+  <Array Name="variables" Type="Text">
+    <Dim>11</Dim>
+    <Stream><Metalink Format="Text" Delimiter=" \n"/>
+z t r11_out r22_out r33_out r12_re_out r12_im_out r13_re_out r13_im_out r23_re_out r23_im_out 
+    </Stream>
+  </Array>
+  <Array Name="data" Type="double">
+    <Dim>201</Dim>
+    <Dim>100</Dim>
+    <Dim>11</Dim>
+    <Stream><Metalink Format="Binary" UnsignedLong="uint32" precision="double" Type="Remote" Encoding="LittleEndian"/>
+Nlevels_no_dopler_with_z_mg1.dat
+    </Stream>
+  </Array>
+</XSIL>
+</simulation>
diff --git a/xmds2/Shahriar_system/Readme b/xmds2/Shahriar_system/Readme
new file mode 100644
index 0000000..c380a31
--- /dev/null
+++ b/xmds2/Shahriar_system/Readme
@@ -0,0 +1,7 @@
+Usual EIT only two fields probe and pump
+./Shahriar_system.run --E2o 30e6 --E3o 0 --d1 0 --da 0 --delta 0
+
+Raman (far detuned)  EIT only two fields probe and pump
+!./Shahriar_system.run --E2o 30e6 --E3o 0 --d1 50e6 --da 50e6 --delta 0
+
+
diff --git a/xmds2/Shahriar_system/Shahriar_system.xmds b/xmds2/Shahriar_system/Shahriar_system.xmds
index 0109168..0d30a32 100644
--- a/xmds2/Shahriar_system/Shahriar_system.xmds
+++ b/xmds2/Shahriar_system/Shahriar_system.xmds
@@ -3,30 +3,29 @@
 
 	<name>Shahriar_system</name>
 
-	<author>Eugeniy Mikhailov</author>
-	<author>Simon Rochester</author>
+	<author>Eugeniy Mikhailov, Simon Rochester</author>
 	<description>
 		License GPL.
 
-		Solving 4 level atom in N-field configuration,
+		Solving 3 level atom in double drive configuration
+		after Shahriar paper about white cavity
 		with field propagation along spatial axis Z
-		no Doppler broadening
+		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).
+		All fields detuned from upper level i.e. Raman configuration
 
-		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>    \
-		*                       \
+		*
+		*         .....
+		*           /   ....
+		*          / .... \
+		*         /   /    \
+		*        /   /-------- |3>  
+		*       E3  /        \
+		*      /   E2         \ 
+		*     /   /            \ E1
+		*    ------ |2>         \
+		*                        \
 		*                    ------- |1>
 		*
 
@@ -55,30 +54,36 @@
 
 				// repopulation rate (atoms flying in/out the laser beam)  in [1/s]
 				const double gt=0.01/2 *(2*M_PI*1e6);
-				// Natural linewidth  of j's level in [1/s]
-				const double G3=2*2.7 *(2*M_PI*1e6);
-				const double G4=2*3.0 *(2*M_PI*1e6);
+				// Natural linewidth  of the upper level [1/s]
+				const double G=2*6 *(2*M_PI*1e6);
+
+				// incoherent pumping rate from level |1> to |3> in [1/s]
+				const double gp=2*.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; // Complex conjugated Rabi frequencies
 
-				complex  r21, r31, r41, r32, r42, r43, r44;  // density matrix elements
+				complex  r21, r31, r32;  // density matrix elements
 			]]>
 		</globals>
 		<benchmark />
     <arguments>
 			<!-- Rabi frequency divided by 2 in [1/s] -->
+			<!--probe-->
       <argument name="E1o" type="real" default_value="0.0025*(2*M_PI*1e6)" />
-      <argument name="E2o" type="real" default_value="6*(2*M_PI*1e6)" />
-      <argument name="E3o" type="real" default_value="0.0*(2*M_PI*1e6)" />
+			<!--pump fields-->
+      <argument name="E2o" type="real" default_value="1.0*(2*M_PI*1e6)" />
+      <argument name="E3o" type="real" default_value="1.0*(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" />
+			<!-- probe field detuning-->
+      <argument name="d1"  type="real" default_value="12*(2*M_PI*1e6)" />
+			<!-- averaged detuning of pump fields i.e. mid point -->
+      <argument name="da"  type="real" default_value="12*(2*M_PI*1e6)" />
+			<!-- detuning of pump fields with respect to each other -->
+      <argument name="delta"  type="real" default_value="6*(2*M_PI*1e6)" />
     </arguments>
 		<bing />
 		<fftw plan="patient" />
@@ -109,13 +114,10 @@
 	</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>
+		<components>r11 r22 r33 r12 r13  r23 </components>
 		<!--
 				 note one of the level population is redundant since
-				 r11+r22+r33+r44=1
-				 so r11 is missing
+				 r11+r22+r33=1
 		-->
 		<initialisation>
 			<![CDATA[
@@ -128,17 +130,28 @@
 			// 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;
+			r11 = 1; r22 = 0; r33 = 0;
+			r12 = 0; r13 = 0;
+			r23 = 0;
 			]]>
 		</initialisation>
 	</vector>
 
+	<vector name="pump_detunings" type="complex" initial_space="t">
+		<components> d dc </components>
+		<!--dc is probably redundant since it is just complex conjugate of d-->
+		<initialisation>
+		  <![CDATA[
+			d  = exp( i*t*delta );
+			dc = conj(d);
+		  ]]>
+		</initialisation>
+	</vector>	
+
+
 	<sequence>
 		<!--For this set of conditions ARK45 is faster than ARK89-->
-		<integrate algorithm="ARK45" tolerance="1e-5" interval="7e-2">
+		<integrate algorithm="ARK45" tolerance="1e-5" interval="4e-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-->
 
@@ -147,13 +160,12 @@
 			<operators>
         <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
 					<integration_vectors>density_matrix</integration_vectors>
-          <dependencies>E_field</dependencies>
+          <dependencies>E_field pump_detunings</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;
+							r11 = 1; r22 = 0; r33 = 0; 
+							r12 = 0; r13 = 0; 
+							r23 = 0; 
             ]]>
           </boundary_condition>
 					<![CDATA[
@@ -163,23 +175,16 @@
 
 						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 + G3*r33*R31 + G4*r44*R41 + i*(-(E1*r13) + E1c*r31);
-						dr12_dt = -(gt*r12) + i*((-delta1 + delta2)*r12 - E2*r13 - E3*r14 + E1c*r32);
-						dr13_dt = -((G3 + 2*gt)*r13)/2. + i*(-(delta1*r13) - E1c*r11 - E2c*r12 + E1c*r33);
-						dr14_dt = -((G4 + 2*gt)*r14)/2. + i*(-((delta1 - delta2 + delta3)*r14) - E3c*r12 + E1c*r34);
-						dr22_dt = gt/2. - gt*r22 + G3*r33*R32 + G4*r44*R42 + i*(-(E2*r23) - E3*r24 + E2c*r32 + E3c*r42);
-						dr23_dt = -((G3 + 2*gt)*r23)/2. + i*(-(delta2*r23) - E1c*r21 - E2c*r22 + E2c*r33 + E3c*r43);
-						dr24_dt = -((G4 + 2*gt)*r24)/2. + i*(-(delta3*r24) - E3c*r22 + E2c*r34 + E3c*r44);
-						dr33_dt = -((G3 + gt)*r33) + i*(E1*r13 + E2*r23 - E1c*r31 - E2c*r32);
-						dr34_dt = -((G3 + G4 + 2*gt)*r34)/2. + i*((delta2 - delta3)*r34 + E1*r14 + E2*r24 - E3c*r32);
-						dr44_dt = -((G4 + gt)*r44) + i*(E3*r24 - E3c*r42);
+						dr11_dt = gt - 2*(gp + gt)*r11 - E1*i*r13 + E1c*i*r31 + G*r33;
+						dr12_dt = (-gp - 2*gt - d1*i + da*i)*r12 - i*(E2*d + E3*dc)*r13 + E1c*i*r32;
+						dr13_dt = -(E1c*i*r11) - i*(E3c*d + E2c*dc)*r12 + (-G - gp - 2*gt - d1*i)*r13 + E1c*i*r33;
+						dr22_dt = gt - 2*gt*r22 - i*(E2*d + E3*dc)*r23 + i*(E3c*d + E2c*dc)*r32 + G*r33;
+						dr23_dt = -(E1c*i*r21) - i*(E3c*d + E2c*dc)*r22 + (-G - 2*gt - da*i)*r23 + i*(E3c*d + E2c*dc)*r33;
+						dr33_dt = 2*gp*r11 + E1*i*r13 + i*(E2*d + E3*dc)*r23 - E1c*i*r31 - i*(E3c*d + E2c*dc)*r32 - 2*(G + gt)*r33;
+
 					]]>
         </operator>
 				<operator kind="ex" constant="yes">
@@ -193,17 +198,14 @@
           <![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] ;
+					dE3_dz = i*eta*conj(r23) -Lt[E3] ;
           ]]>
 			</operators>
 		</integrate>
 	</sequence>
 
-
-
-
 	<!-- The output to generate -->
-	<output format="binary" filename="Nlevels_no_dopler_with_z.xsil">
+	<output format="binary" filename="Shahriar_system.xsil">
 		<group>
       <sampling basis="t(1000)" initial_sample="yes">
 				<dependencies>E_field</dependencies>
@@ -220,30 +222,22 @@
       <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
+					r11_out r22_out r33_out
+					r12_re_out r12_im_out r13_re_out r13_im_out
+					                      r23_re_out r23_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>
diff --git a/xmds2/Shahriar_system/pp.m b/xmds2/Shahriar_system/pp.m
new file mode 100644
index 0000000..fcb7e2a
--- /dev/null
+++ b/xmds2/Shahriar_system/pp.m
@@ -0,0 +1,106 @@
+Shahriar_system
+
+%% field propagation
+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')
+ylabel('t (uS)')
+zlabel('I_1')
+title('I_1')
+subplot(1,3,2); imagesc(z_1, t_1, I2_out_1); colorbar
+xlabel('z')
+ylabel('t (uS)')
+zlabel('I_2')
+title('I_2')
+subplot(1,3,3); imagesc(z_1, t_1, I3_out_1); colorbar
+xlabel('z')
+ylabel('t (uS)')
+zlabel('I_3')
+title('I_3')
+
+
+
+
+%% fields before and after the cell
+figure(2)
+subplot(1,3,1);
+plot( ...
+	t_1,I1_out_1(:,1),'-;before;',  ...
+	t_1,I1_out_1(:,end), '-;after;' ...
+	)
+xlabel('t')
+ylabel('I_1')
+title('I_1 propagation')
+subplot(1,3,2);
+plot( ...
+	t_1,I2_out_1(:,1),'-;before;',  ...
+	t_1,I2_out_1(:,end), '-;after;' ...
+	)
+xlabel('t')
+ylabel('I_2')
+title('I_2 propagation')
+subplot(1,3,3);
+plot( ...
+	t_1,I3_out_1(:,1),'-;before;',  ...
+	t_1,I3_out_1(:,end), '-;after;' ...
+	)
+xlabel('t')
+ylabel('I_3')
+title('I_3 propagation')
+
+%return;
+
+%% all density matrix elements in one plot
+% diagonal populations, 
+% upper triangle real part of coherences, 
+% lower diagonal imaginary part of coherences
+figure(3)
+subplot(3,3,1); imagesc (z_2, t_2, r11_out_2);  caxis([0,1]); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('rho_{11}')
+title('rho_{11}')
+subplot(3,3,5); imagesc (z_2, t_2, r22_out_2);  caxis([0,1]); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('rho_{22}')
+title('rho_{22}')
+subplot(3,3,9); imagesc (z_2, t_2, r33_out_2);  caxis([0,1]); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('rho_{33}')
+title('rho_{33}')
+% real parts of coherences
+subplot(3,3,2); imagesc(z_2, t_2, r12_re_out_2); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('Real(rho_{12})')
+title('Real(rho_{12})')
+subplot(3,3,3); imagesc(z_2, t_2, r13_re_out_2); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('Real(rho_{13})')
+title('Real(rho_{13})')
+subplot(3,3,6); imagesc(z_2, t_2, r23_re_out_2); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('Real(rho_{23})')
+title('Real(rho_{23})')
+% imaginary parts of coherences
+subplot(3,3,4); imagesc(z_2, t_2, r12_im_out_2); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('Imag(rho_{12})')
+title('Imag(rho_{12})')
+subplot(3,3,7); imagesc(z_2, t_2, r13_im_out_2); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('Imag(rho_{13})')
+title('Imag(rho_{13})')
+subplot(3,3,8); imagesc(z_2, t_2, r23_im_out_2); colorbar
+xlabel('z')
+ylabel('t')
+zlabel('Imag(rho_{23})')
+title('Imag(rho_{23})')
+