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