diff options
6 files changed, 666 insertions, 0 deletions
diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/Makefile b/xmds2/Nlevels_with_doppler_with_z_4wm/Makefile new file mode 100644 index 0000000..578f2b0 --- /dev/null +++ b/xmds2/Nlevels_with_doppler_with_z_4wm/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_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds new file mode 100644 index 0000000..8faa04d --- /dev/null +++ b/xmds2/Nlevels_with_doppler_with_z_4wm/Nlevels_with_doppler_with_z_4wm.xmds @@ -0,0 +1,322 @@ +<?xml version="1.0"?> +<simulation xmds-version="2"> + + <name>Nlevels_with_doppler_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', 't', and 'v' to have dimensions [m], [s], and [m/s] --> + <geometry> + <propagation_dimension> z </propagation_dimension> + <transverse_dimensions> + <dimension name="t" lattice="1000" domain="(-14.0e-7, 4.0e-7)" /> + <dimension name="v" lattice="1" domain="(-1000, 1000)" /> + </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_with_doppler_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_with_doppler_with_z_4wm/Readme b/xmds2/Nlevels_with_doppler_with_z_4wm/Readme new file mode 100644 index 0000000..aa04999 --- /dev/null +++ b/xmds2/Nlevels_with_doppler_with_z_4wm/Readme @@ -0,0 +1,14 @@ +Fast light of field E2 +./Nlevels_no_dopler_with_z_4wm.run --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=3.8e7 --E4o=6.3e4 + +Slow light +./Nlevels_no_dopler_with_z_4wm.run --delta1=0 --delta2=0 --delta3=0 --E1o=1.9e7 --E2o=3.1e5 --E3o=0 --E4o=0 + +Sigar shaped propagation +./Nlevels_no_dopler_with_z_4wm.run --delta1=0 --delta2=0 --delta3=0 --E1o=6e6 --E2o=3e5 --E3o=6.0e6 --E4o=3e5 + +Fast light to Slow light switch +./Nlevels_no_dopler_with_z_4wm.run --delta1=0 --delta2=0 --delta3=0 --E1o=2e7 --E2o=3e3 --E3o=6e6 --E4o=3e3 + + + diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/plot_fields_propagation_I2.gp b/xmds2/Nlevels_with_doppler_with_z_4wm/plot_fields_propagation_I2.gp new file mode 100644 index 0000000..b6721d4 --- /dev/null +++ b/xmds2/Nlevels_with_doppler_with_z_4wm/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.4:0.4] 'I2.dat' diff --git a/xmds2/Nlevels_with_doppler_with_z_4wm/pp.m b/xmds2/Nlevels_with_doppler_with_z_4wm/pp.m new file mode 100644 index 0000000..b88fa6e --- /dev/null +++ b/xmds2/Nlevels_with_doppler_with_z_4wm/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_with_doppler_with_z_4wm/pp_I2.m b/xmds2/Nlevels_with_doppler_with_z_4wm/pp_I2.m new file mode 100644 index 0000000..45d4913 --- /dev/null +++ b/xmds2/Nlevels_with_doppler_with_z_4wm/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') + |