Nlevels_with_MOR

%% output E fields to intensities
IdL_out_1 = abs(EdLaR_1 +1i*EdLaI_1).^2;
IdR_out_1 = abs(EdRaR_1 +1i*EdRaI_1).^2;
IpL_out_1 = abs(EpLaR_1 +1i*EpLaI_1).^2;
IpR_out_1 = abs(EpRaR_1 +1i*EpRaI_1).^2;

EpR = EpRaR_1 + 1i*EpRaI_1;
EpL = EpLaR_1 + 1i*EpLaI_1;

Ieit = abs(EpR).^2 + abs(EpL).^2;
Epos = 1/sqrt(2)*(EpR + 1i*EpL);
Eneg = 1/sqrt(2)*(EpR - 1i*EpL);
Imor = abs(Epos).^2 - abs(Eneg).^2;

%% Maps are not so useful with too many pixels in the image
% Here we limit number of image points
Npx = 800;
Npy=600;

Ndx = length(z_1);
Ndy = length(t_1);

skip_xp = max(1, floor(Ndx/Npx) );
skip_yp = max(1, floor(Ndy/Npy) );
ind_x = 1:skip_xp:Ndx;
ind_y = 1:skip_yp:Ndy;


%% 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(ind_x), t_1(ind_y), IdL_out_1(ind_x, ind_y).'); colorbar
xlabel('z (cm)')
ylabel('t (uS)')
zlabel('I_{dL}')
title('I_{dL}')
subplot(2,2,3); imagesc(z_1(ind_x), t_1(ind_y), IpL_out_1(ind_x, ind_y).'); colorbar
xlabel('z (cm)')
ylabel('t (uS)')
zlabel('I_{pL}')
title('I_{pL}')
subplot(2,2,2); imagesc(z_1(ind_x), t_1(ind_y), IdR_out_1(ind_x, ind_y).'); colorbar
xlabel('z (cm)')
ylabel('t (uS)')
zlabel('I_{dR}')
title('I_{dR}')
subplot(2,2,4); imagesc(z_1(ind_x), t_1(ind_y), IpR_out_1(ind_x, ind_y).'); 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,3);
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,2);
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')

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')
legend('before', 'after')

%% EIT and MOR fields before and after the cell
figure(6)
subplot(2,2,1);
plot( ...
	t_1, (Ieit(1,:)'),   ...
	t_1, (Ieit(end,:)'),'LineWidth', 4 ...
	)
xlabel('t (uS)')
ylabel('I_{eit} (1/s)^2')
title('I_{eit} before and after cell')
legend('before', 'after')

%%
subplot(2,2,3);
plot( ...
	t_1, (Imor(1,:)'), ...
	t_1, (Imor(end,:)'), 'linewidth', 4 ...
	)
xlabel('t (uS)')
ylabel('I_{mor} (1/s)^2')
title('I_{mor} before and after cell')
legend('before', 'after')

%%
subplot(2,2,2);
imagesc(z_1(ind_x), t_1(ind_y), Ieit(ind_x, ind_y).'); colorbar
xlabel('z (cm)')
ylabel('t (uS)')
zlabel('I_{eit}')
title('I_{eit}')

subplot(2,2,4);
imagesc(z_1(ind_x), t_1(ind_y), Imor(ind_x, ind_y).'); colorbar
xlabel('z (cm)')
ylabel('t (uS)')
zlabel('I_{mor}')

title('I_{mor}')

%% EIT and MOR signals after the cell side by side
figure(3)
plot( ...
	t_1, (Ieit(end,:)'), ...
	t_1, (Imor(end,:)')
	)
xlabel('t (uS)')
ylabel('I_{eit} or I_{MOR} (1/s)^2')
title('I_{eit} and I_{mor} after cell')
legend('EIT', 'MOR')


%%
return

[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')


%%
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),flip(IpR_in_norm(indx)),  ...
	t_1(indx),flip(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})')