full simulation of 5 level MOR to separate folder

3 files changed, 723 insertions, 0 deletions

diff --git a/MOR_5_levels_with_doppler_and_propagation/Makefile b/MOR_5_levels_with_doppler_and_propagation/Makefile
new file mode 100644
index 0000000..ade8cb0
--- /dev/null
+++ b/MOR_5_levels_with_doppler_and_propagation/Makefile
@@ -0,0 +1,55 @@
+###  -*- 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 = xsil2graphics2
+
+all: $(RUN_FILES) 
+
+%.run: %.xmds
+	$(XMDS) $<
+	mv $(patsubst %.xmds,%,$<) $@
+
+%.xsil: %.run
+	./$<  --E1o=0 --E3o=0
+
+%.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)
+	rm -f $(eps_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/MOR_5_levels_with_doppler_and_propagation/Nlevels_with_MOR.xmds b/MOR_5_levels_with_doppler_and_propagation/Nlevels_with_MOR.xmds
new file mode 100644
index 0000000..99d88a0
--- /dev/null
+++ b/MOR_5_levels_with_doppler_and_propagation/Nlevels_with_MOR.xmds
@@ -0,0 +1,481 @@
+<?xml version="1.0"?>
+<simulation xmds-version="2">
+
+	<name>Nlevels_with_MOR</name>
+
+	<author>Eugeniy Mikhailov, M. Guidry</author>
+	<description>
+		License GPL.
+
+		Solving split 3-level atom
+		with field propagation along spatial axis Z
+		with Doppler broadening.
+
+
+		*
+		*                  -------- |a>  
+		*                 /  / \  \
+		*           EdL  /  /   \  \ 
+		*               /  /     \  \   
+		*   |c> ----------/---    \  \ EpR
+		*                / EpL     \  \
+		*   |C> --------------      \  \
+		*                       EdR  \  \
+		*                             \  \
+		*                        ------\------ |b>
+		*                               \
+		*                        ------------- |B>
+		*
+		*
+
+
+		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
+		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 k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K]
+				const double lambda=794.7e-9; //wavelength in m
+				const double Kvec = 2*M_PI/lambda;  // k-vector
+				const double Gamma_super=6*(2*M_PI*1e6);  // characteristic decay rate of upper level used for eta  calculations expressed in [1/s]
+				// eta will be calculated in the <arguments> section
+				double eta = 0;  // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				
+				//  ---------  Atom and cell properties -------------------------
+				// range of Maxwell distribution atomic velocities
+				const double mass = (86.909180527 * 1.660538921e-27); // atom mass in [kg] 
+				// above mass expression is written as (expression is isotopic_mass * atomic_mass_unit)
+
+				// Average sqrt(v^2) in Maxwell distribution for one dimension
+				// Maxwell related parameters will be calculated in <arguments> section
+				double v_thermal_averaged=0;
+				// Maxwell distribution velocities range to take in account in [m/s]
+				double V_maxwell_min = 0, V_maxwell_max = 0;
+
+				// 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 Ga=3.0 *(2*M_PI*1e6);
+				//const double G4=3.0 *(2*M_PI*1e6);
+
+				complex g = Gamma_super;
+				complex gbc = (2*M_PI)*1e3; // 1/s untits
+				const complex Split = 0;
+				const complex noise = 0;
+				
+				complex Gab, GAB, Gca, GCA, Gcb, GCB;
+
+				// total decay of i-th level branching ratios. Rij branching of i-th level to j-th
+				//const double Rab=0.5, Rac=0.5;
+
+
+				complex EdLac, EdRac, EpLac, EpRac; 
+
+				// inner use variables 
+				double probability_v; // will be used as p(v) in Maxwell distribution
+
+			]]>
+		</globals>
+		<validation kind="run-time"/> <!--allows to put ranges as variables-->
+		<benchmark />
+    <arguments>
+			<!-- Rabi frequency divided by 2 in [1/s] -->
+      <argument name="EdLo" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
+      <argument name="EdRo" type="real" default_value="2*1.5*(2*M_PI*1e6)" />
+
+      <argument name="EpLo" type="real" default_value="0.05*(2*M_PI*1e6)" /> 
+      <argument name="EpRo" type="real" default_value="0.05*(2*M_PI*1e6)" />
+ 
+
+			<!-- Fields detuning in [1/s] -->
+      <argument name="delta_dL"  type="real" default_value="0.0" />
+      <argument name="delta_pL"  type="real" default_value="0.0" />
+      <argument name="delta_dR"  type="real" default_value="0.0" />
+      <argument name="delta_pR"  type="real" default_value="0.0" />
+			<!--Pulse duration/width [s] -->
+			<argument name="Pwidth"  type="real" default_value="0.1e-6" />
+			<!--  Atom and cell properties -->
+			<!--Cell length [m] -->
+			<argument name="Lcell"  type="real" default_value="1.5e-2" />
+			<!--Density of atoms [1/m^3] -->
+			<argument name="Ndens"  type="real" default_value="1e15" />
+			<!--Atoms temperature [K] -->
+			<!--TODO: looks like Temperature > 10 K knocks solver, 
+					 I am guessing detunings are too large and thus it became a stiff equation-->
+			<!--! make sure it is not equal to zero!-->
+			<argument name="Temperature"  type="real" default_value="5" />
+			<!-- This will be executed after arguments/parameters are parsed -->
+			<!-- Read the code Luke: took me a while of reading the xmds2 sources to find it -->
+			<![CDATA[
+				// Average sqrt(v^2) in Maxwell distribution for one dimension
+				if (Temperature == 0)
+					_LOG(_ERROR_LOG_LEVEL, "ERROR: Temperature should be >0 to provide range for Maxwell velocity distribution\n");
+				v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); 
+				// Maxwell distribution velocities range to take in account in [m/s]
+				// there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0)
+				V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; 
+
+				// eta constant in the wave equation for Rabi frequency. Units are [1/(m s)]
+				eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI;
+			]]>
+    </arguments>
+		<bing />
+		<diagnostics /> 
+		<fftw plan="estimate" threads="1" />
+		<!-- I don't see any speed up on 6 core CPU even if use threads="6" -->
+		<openmp />
+		<auto_vectorise />
+		<halt_non_finite />
+	</features>
+
+	<!-- 'z', 't', and 'v'  to have dimensions [m], [s], and [m/s]   -->
+	<geometry>
+		<propagation_dimension> z </propagation_dimension>
+		<transverse_dimensions>
+			<!-- IMPORTANT: looks like having a lot of points in time helps with convergence.
+					 I suspect that time spacing should be small enough to catch
+					 all pulse harmonics and more importantly 1/dt should be larger than
+					 the largest detuning (including Doppler shifts).
+					 Unfortunately calculation time is proportional to lattice size
+					 so we cannot just blindly increase it.
+					 Some rules of thumb:  
+						* lattice="1000"   domain="(-1e-6, 1e-6)" 
+							was good enough detunings up to 155 MHz (980 rad/s) notice that 1/dt=500 MHz
+						* lattice="10000"   domain="(-1e-6, 1e-6)" 
+							works for Doppler averaging in up to 400K for Rb when lasers are zero detuned
+			 -->
+			<dimension name="t"   lattice="10000"   domain="(-1e-6, 1e-6)" />
+			<dimension name="v"   lattice="2"   domain="(V_maxwell_min, V_maxwell_max)" />
+		</transverse_dimensions>
+	</geometry>
+
+	<!-- Rabi frequency --> 
+	<vector name="E_field" type="complex" initial_space="t">
+		<components>EdL EdR EpL EpR</components>
+		<initialisation>
+			<![CDATA[
+			// Initial (at starting 'z' position) electromagnetic field does not depend on detuning
+			// as well as time
+			EdL=EdLo*exp(-pow( ((t-0.0)/Pwidth),2) );
+			EdR=EdRo*exp(-pow( ((t-0.0)/Pwidth),2) );
+			EpL=EpLo*exp(-pow( ((t-0.0)/Pwidth),2) );
+			EpR=EpRo*exp(-pow( ((t-0.0)/Pwidth),2) );
+
+			//EpL = EpLo*(1+0.01*2*(((double)rand() / (double)RAND_MAX)-0.5));	
+			//EpR = EpRo*(1+0.01*2*(((double)rand() / (double)RAND_MAX)-0.5));
+
+			//E2=E2o*exp(-pow( ((t-0.0)/Pwidth),2) );
+			]]>
+		</initialisation>
+	</vector>
+
+	<!--Maxwell distribution probability p(v)-->
+	<computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real">
+		<components>probability_v</components>
+		<evaluation>
+			<![CDATA[
+			// TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during  the script lifetime since 'v' is fixed.
+			probability_v=1.0/(v_thermal_averaged*sqrt(2*M_PI)) * exp( - mod2(v/v_thermal_averaged)/2.0 ); 
+			]]>
+		</evaluation>
+	</computed_vector>
+
+	<!--Maxwell distribution norm sum(p(v))
+			 Needed since we sum over the grid instead of true integral,
+			 we also have finite cut off velocities-->
+	<computed_vector name="Maxwell_distribution_probabilities_norm" dimensions="" type="real">
+		<components>probability_v_norm</components>
+		<evaluation>
+			<dependencies basis="v">Maxwell_distribution_probabilities</dependencies>
+			<![CDATA[
+			// TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during  the script lifetime since 'v' is fixed.
+			probability_v_norm=probability_v;
+			]]>
+		</evaluation>
+	</computed_vector>
+
+
+	<!-- Averaged across Maxwell distribution fields amplitudes -->
+	<computed_vector name="E_field_avgd" dimensions="t" type="complex">
+		<components>EdLa EdRa EpLa EpRa</components>
+		<evaluation>
+			<dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
+			<![CDATA[
+			double prob_v_normalized=probability_v/probability_v_norm;
+
+			EdLa=EdL*prob_v_normalized;
+			EdRa=EdR*prob_v_normalized;
+			EpLa=EpL*prob_v_normalized;
+			EpRa=EpR*prob_v_normalized;
+			]]>
+		</evaluation>
+	</computed_vector>
+
+	<!-- Averaged across Maxwell distribution density matrix components -->
+	<computed_vector name="density_matrix_averaged" dimensions="t" type="complex">
+		<components>rbb_av rBB_av rcc_av rCC_av raa_av rcb_av rab_av rca_av rCB_av rAB_av rCA_av</components>
+		<evaluation>
+			<dependencies basis="v">density_matrix Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies>
+			<![CDATA[
+			double prob_v_normalized=probability_v/probability_v_norm;
+
+			rbb_av=rbb*prob_v_normalized;
+			rBB_av=rBB*prob_v_normalized;
+			rcc_av=rcc*prob_v_normalized;
+			rCC_av=rCC*prob_v_normalized;
+			raa_av=raa*prob_v_normalized;
+
+			rcb_av=rcb*prob_v_normalized;
+			rab_av=rab*prob_v_normalized;
+			rca_av=rca*prob_v_normalized;
+
+			rCB_av = rCB*prob_v_normalized;
+			rAB_av = rAB*prob_v_normalized;
+			rCA_av = rCA*prob_v_normalized;
+
+			]]>
+		</evaluation>
+	</computed_vector>
+
+
+	<vector name="density_matrix" type="complex" initial_space="t">
+		<components>rbb rBB rcc rCC raa rcb rab rca rCB rAB rCA</components>
+		<initialisation>
+			<!--This sets boundary condition at all times and left border of z (i.e. z=0)-->
+				<dependencies>E_field_avgd</dependencies>
+				<![CDATA[
+						EdLac = conj(EdLa);
+						EdRac = conj(EdRa);
+						EpLac = conj(EpLa);
+						EpRac = conj(EpRa);
+
+				rbb = 0.25; rcc = 0.25; rBB = 0.25; rCC = 0.25;
+				raa = 0;
+				rcb = 0; rab = 0; rca = 0;
+				rCB = 0; rAB = 0; rCA = 0;
+
+				]]>
+		</initialisation>
+	</vector>
+
+	<sequence>
+		<!--For this set of conditions ARK45 is faster than ARK89-->
+		<!--ARK45 is good for small detuning when all frequency like term are close to zero-->
+		<integrate algorithm="ARK45" tolerance="1e-5" interval="Lcell"> 
+		<!--<integrate algorithm="SI" steps="200" interval="Lcell"> -->
+		<!--RK4 is good for large detunings when frequency like term are big, it does not try to be too smart about adaptive step which ARK seems to make too small-->
+		<!--When ARK45 works it about 3 times faster then RK4 with 1000 steps-->
+		<!--<integrate algorithm="RK4" steps="100" 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>100 100</samples>
+			<!--Use the next line for debuging to see velocity dependence. Uncomment/switch on output groups 3,4-->
+			<!--<samples>100 100 100 100</samples>--> 
+			<operators>
+        <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t">
+					<integration_vectors>density_matrix</integration_vectors>
+          <dependencies>E_field_avgd</dependencies>
+          <boundary_condition kind="left">
+						<!--This set boundary condition at all 'z'  and left border of 't' (i.e. min(t))-->
+          </boundary_condition>
+					<![CDATA[
+
+					Gab=g+i*(delta_dL+delta_pR+0*noise);
+					GAB=g+i*(delta_pL+delta_dR+0*noise);
+
+					Gca=g-i*(delta_dL);
+					GCA=g-i*(delta_pL);
+
+					Gcb=gbc+i*(Split + delta_pR+0*noise);
+					GCB=gbc+i*(-Split + delta_pL+0*noise);
+
+					complex rba=conj(rab);
+					complex rac=conj(rca);
+					complex rbc=conj(rcb);
+					complex rBA=conj(rAB);
+					complex rAC=conj(rCA);
+					complex rBC=conj(rCB);
+
+					draa_dt = -i*EpRac*rab+i*EpRa*rba-i*EdLac*rac+i*EdLa*rca-2*g*raa
+					          -i*EdRac*rAB+i*EdRa*rBA-i*EpLac*rAC+i*EpLa*rCA;
+
+					drbb_dt =  i*EpRac*rab-i*EpRa*rba+g*raa-gbc*rbb+gbc*rcc;
+					drBB_dt =  i*EdRac*rAB-i*EdRa*rBA+g*raa-gbc*rBB+gbc*rCC;
+
+					drcc_dt =  i*EdLac*rac-i*EdLa*rca+g*raa-gbc*rcc+gbc*rbb;
+					drCC_dt =  i*EpLac*rAC-i*EpLa*rCA+g*raa-gbc*rCC+gbc*rBB;
+
+					drab_dt = -Gab*rab+i*EpRa*(rbb-raa)+i*EdLa*rcb;
+					drca_dt = -Gca*rca+i*EdLac*(raa-rcc)-i*EpRac*rcb;
+					drcb_dt = -Gcb*rcb-i*EpRa*rca+i*EdLac*rab;
+
+					drAB_dt = -Gab*rAB+i*EdRa*(rBB-raa)+i*EpLa*rCB;
+					drCA_dt = -Gca*rCA+i*EpLac*(raa-rCC)-i*EdRac*rCB;
+					drCB_dt = -GCB*rCB-i*EdRa*rCA+i*EpLac*rAB;
+
+
+					]]>
+        </operator>
+				<!--
+							 According to xmds2 docs operator kind="ip" should be faster
+							 but our codes runs about 5% to 10% slower with it.
+							 Maybe because we very close to the stiff condition so I use "ex" kind
+							 <operator kind="ip" constant="yes">
+					 -->
+				<operator kind="ex" constant="yes" type="imaginary">
+					<operator_names>Lt</operator_names>
+					<![CDATA[
+					Lt = -i/c*kt;
+					]]>
+				</operator>
+        <integration_vectors>E_field</integration_vectors>
+				<dependencies>density_matrix</dependencies>
+          <![CDATA[
+					dEdL_dz = i*eta*(rca) +0*Lt[EdL] ;
+					dEdR_dz = i*eta*(rAB) +0*Lt[EdR] ;
+					dEpL_dz = i*eta*(rCA) +0*Lt[EpL] ;
+					dEpR_dz = i*eta*(rab) +0*Lt[EpR] ;
+ 
+          ]]>
+			</operators>
+		</integrate>
+	</sequence>
+
+
+
+	<!-- The output to generate -->
+	<output format="binary" filename="Nlevels_with_MOR.xsil">
+		<group>
+      <sampling basis="t(1000) " initial_sample="yes">
+				<dependencies>E_field_avgd</dependencies>
+				<moments>IdL_out IpL_out IdR_out IpR_out</moments>
+				<![CDATA[
+				IdL_out = mod2(EdLa);
+				IpL_out = mod2(EpLa);
+				IdR_out = mod2(EdRa);
+				IpR_out = mod2(EpRa);
+				]]>
+			</sampling>
+		</group>
+
+		<group>
+      <sampling basis="t(100) v(2)" initial_sample="yes">
+				<dependencies>density_matrix_averaged</dependencies>
+				<moments>
+					rbb_out rBB_out
+					rcc_out rCC_out
+					raa_out 
+					rcb_re_out rcb_im_out
+ 					rCB_re_out rCB_im_out
+					rab_re_out rab_im_out 
+					rAB_re_out rAB_im_out
+					rca_re_out rca_im_out 
+					rCA_re_out rCA_im_out
+				</moments>
+				<![CDATA[
+				// populations output 
+				rbb_out = rbb_av.Re();
+				rBB_out = rBB_av.Re();
+
+				rcc_out = rcc_av.Re();	
+				rCC_out = rCC_av.Re();
+			
+				raa_out = raa_av.Re();
+
+				// coherences output 
+				rcb_re_out = rcb_av.Re();
+				rcb_im_out = rcb_av.Im();
+
+				rCB_re_out = rCB_av.Re();
+				rCB_im_out = rCB_av.Im();
+				
+				rab_re_out = rab_av.Re();
+				rab_im_out = rab_av.Im();
+
+				rAB_re_out = rAB_av.Re();
+				rAB_im_out = rAB_av.Im();
+
+				rca_re_out = rca_av.Re();
+				rca_im_out = rca_av.Im();
+
+				rCA_re_out = rCA_av.Re();
+				rCA_im_out = rCA_av.Im();
+
+				]]>
+			</sampling>
+		</group>
+
+		<!-- use the following two groups only for debuging 
+				 otherwise they are quite useless and have to much information 
+				 in 3D space (z,t,v) -->
+		<!--
+		<group>
+      <sampling basis="t(100) v(10)" initial_sample="yes">
+				<dependencies>E_field</dependencies>
+				<moments>I1_out_v I2_out_v I3_out_v I4_out_v</moments>
+				<![CDATA[
+				// light field intensity distribution in velocity subgroups 
+				I1_out_v = mod2(E1);
+				I2_out_v = mod2(E2);
+				I3_out_v = mod2(E3);
+				I4_out_v = mod2(E4);
+				]]>
+			</sampling>
+		</group>
+
+		<group>
+      <sampling basis="t(100) v(10)" initial_sample="yes">
+				<dependencies>density_matrix</dependencies>
+				<moments>
+					rbb_out_v rcc_out_v raa_out_v r44_out_v 
+					rbc_re_out_v rbc_im_out_v rba_re_out_v rba_im_out_v r14_re_out_v r14_im_out_v
+					                      rca_re_out_v rca_im_out_v r24_re_out_v r24_im_out_v 
+					                                            r34_re_out_v r34_im_out_v
+				</moments>
+				<![CDATA[
+				// density matrix distribution in velocity subgroups 
+				// populations output 
+				rbb_out_v = rbb.Re();
+				rcc_out_v = rcc.Re();
+				raa_out_v = raa.Re();
+				r44_out_v = r44.Re();
+				// coherences output 
+				rbc_re_out_v = rbc.Re();
+				rbc_im_out_v = rbc.Im();
+				rba_re_out_v = rba.Re();
+				rba_im_out_v = rba.Im();
+				r14_re_out_v = r14.Re();
+				r14_im_out_v = r14.Im();
+				rca_re_out_v = rca.Re();
+				rca_im_out_v = rca.Im();
+				r24_re_out_v = r24.Re();
+				r24_im_out_v = r24.Im();
+				r34_re_out_v = r34.Re();
+				r34_im_out_v = r34.Im();
+				]]>
+			</sampling>
+		</group>
+		-->
+
+	</output>
+
+</simulation>
+	
+<!-- 
+vim: ts=2 sw=2 foldmethod=indent: 
+-->
diff --git a/MOR_5_levels_with_doppler_and_propagation/pp_Nlevels.m b/MOR_5_levels_with_doppler_and_propagation/pp_Nlevels.m
new file mode 100644
index 0000000..0991c97
--- /dev/null
+++ b/MOR_5_levels_with_doppler_and_propagation/pp_Nlevels.m
@@ -0,0 +1,187 @@
+Nlevels_with_MOR
+
+%% 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, IdL_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_{dL}')
+title('I_{dL}')
+subplot(2,2,2); imagesc(z_1, t_1, IpL_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_{pL}')
+title('I_{pL}')
+subplot(2,2,3); imagesc(z_1, t_1, IdR_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_{dR}')
+title('I_{dR}')
+subplot(2,2,4); imagesc(z_1, t_1, IpR_out_1); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('I_{pR}')
+title('I_{pR}')
+
+
+%print('-color','fields_propagation.eps')
+
+
+
+%% fields before and after the cell
+figure(5)
+subplot(2,2,1);
+plot( ...
+	t_1,IdL_out_1(1,:)',   ...
+	t_1,IdL_out_1(end,:)','LineWidth', 4 ...
+	)
+xlabel('t (uS)')
+ylabel('I_{dL} (1/s)^2')
+title('I_{dL} before and after cell')
+legend('before', 'after')
+
+%%
+subplot(2,2,2);
+plot( ...
+	t_1,IpL_out_1(1,:)', ...
+	t_1,IpL_out_1(end,:)', 'linewidth', 4 ...
+	)
+xlabel('t (uS)')
+ylabel('I_{pL} (1/s)^2')
+title('I_{pL} before and after cell')
+legend('before', 'after')
+
+%%
+subplot(2,2,3);
+plot( ...
+	t_1,IdR_out_1(1,:)',  ...
+	t_1,IdR_out_1(end,:)', 'linewidth', 4 ...
+	)
+xlabel('t (uS)')
+ylabel('I_{dR} (1/s)^2')
+title('I_{dR} before and after cell')
+legend('before', 'after')
+
+%%
+
+[b, a]=butter(3, 0.05);
+
+IpR_out_after=IpR_out_1(end,:);
+
+IpR_out_after_filtered=filtfilt(b,a,IpR_out_after);
+settling_time=0.8; %uS
+t_good_indx=t_1> min(t_1 + settling_time);
+
+[m, max_pos_before]=max(IpR_out_1(1,t_good_indx) );
+
+
+[m, max_pos_after]=max(IpR_out_after_filtered(1, t_good_indx));
+
+
+delay_time=t_1(max_pos_after)-t_1(max_pos_before);
+
+
+display( strcat('Second field delay time  = ', num2str(delay_time), ' uS/n'));
+
+%%
+
+%print('-color','fields_before_after_cell.eps')
+
+subplot(2,2,4);
+plot( ...
+	t_1,IpR_out_1(1,:)', ...
+	t_1,IpR_out_1(end,:)', 'linewidth', 4 ...
+	)
+xlabel('t (uS)')
+ylabel('I_{pR} (1/s)^2')
+title('I_{pR} before and after cell')
+
+
+%%
+figure(4)
+IpR_max_in=max(IpR_out_1(1,t_good_indx));
+IpR_max_out=max(IpR_out_1(end, t_good_indx));
+IpR_in_norm=(IpR_out_1(1,:))/IpR_max_in;
+IpR_out_norm=(IpR_out_1(end,:))/IpR_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),IpR_in_norm(indx),  ...
+	t_1(indx),IpR_out_norm(indx), 'linewidth', 4 ...
+	)
+xlim([tmin,tmax]);	
+xlabel('t (uS)')
+ylabel('I_{pR}')
+title('I_{pR} before and after cell normalized')
+%print('-color','probe_before_after_cell.eps')
+legend('before', 'after')
+
+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_1, t_2, rbb_out_2(:,:,1));  caxis([0,1]); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('rho_{bb}')
+title('rho_{bb}')
+
+%%
+subplot(4,4,6); imagesc (z_2, t_2, rcc_out_2(:,:,1));  caxis([0,1]); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('rho_{cc}')
+title('rho_{cc}')
+subplot(4,4,11); imagesc (z_2, t_2, raa_out_2);  caxis([0,1]); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('rho_{aa}')
+title('rho_{aa}')
+
+% real parts of coherences
+subplot(4,4,2); imagesc(z_2, t_2, rcb_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{cb})')
+title('Real(rho_{cb})')
+subplot(4,4,3); imagesc(z_2, t_2, rab_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{ab})')
+title('Real(rho_{ab})')
+
+subplot(4,4,7); imagesc(z_2, t_2, rca_re_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Real(rho_{ca})')
+title('Real(rho_{ca})')
+
+% imaginary parts of coherences
+subplot(4,4,5); imagesc(z_2, t_2, rcb_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{cb})')
+title('Imag(rho_{cb})')
+subplot(4,4,9); imagesc(z_2, t_2, rab_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{ab})')
+title('Imag(rho_{ab})')
+subplot(4,4,10); imagesc(z_2, t_2, rca_im_out_2); colorbar
+xlabel('z (cm)')
+ylabel('t (uS)')
+zlabel('Imag(rho_{ca})')
+title('Imag(rho_{ca})')
+