diff options
Diffstat (limited to 'MOR_5_levels_with_doppler_and_propagation')
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})') + |