% Spatial plot

%% The first section is only for internal heat excitation. 
% For other cases, load data directly and go to the next section
nsteps = 1250;

Data = load('pres2.5mm_below_0.9active_2.0r1G_1250_200.mat');
temp = 1;
ubm1temp(:,:,temp) = Data.ubm1;
uowtemp(temp,:) = Data.uow;
urwtemp(temp,:) = Data.urw;
Data = load('pres2.5mm_above_0.9active_2.0r1G_1250_200.mat');
temp = 2;
ubm1temp(:,:,temp) = Data.ubm1;
uowtemp(temp,:) = Data.uow;
urwtemp(temp,:) = Data.urw;

ubm1 = zeros(441,nsteps);
uow = zeros(1,nsteps);
urw = zeros(1,nsteps);

weight = [0.52 0.48];
for count = 1:length(weight)
    ubm1 = ubm1 + ubm1temp(:,:,count)*weight(count);
    uow  = uow  + uowtemp(count,:)*weight(count);
    urw  = urw  + urwtemp(count,:)*weight(count);
end

%%

nstep = 1250;
lengths = 2*nstep;              % vector length used for ifft
df = 0.2;                   % unit: kHz
freq = linspace(df,nsteps*df,nsteps);

dt = 1/lengths/df;           % unit: msec
t = linspace(0,1/df,lengths+1);

% in cgs units
beta = 385e-6; % 1/K, at 40 C
Cp = 4.178e7; % cm^2/K/s^2, at 37 C
H0 = 0.5e7; % g cm^2/s^3
rhof = 1; % g/cm^3
alpha = 0.3;
mu = 0.00; % fluid viscosity in cgs
c = 1.5d5; % speed of sound in cgs

vbm1 = zeros(size(ubm1));
adjBM = zeros(size(ubm1));
adjOW = zeros(size(uow));

for ic = 1:nsteps
    omega = freq(ic)*1000*2*pi;
    alpha_f = 1 - 1j*omega*4*mu/3/rhof/c^2;
    vbm1(:,ic) = ubm1(:,ic)*(-1i)*omega;
    adjBM(:,ic) = ubm1(:,ic)/rhof/omega^2*alpha_f*alpha*beta*H0/Cp*(1-exp(1i*omega*0.03*1d-3))*(-1i)*omega;
    adjOW(:,ic) = uow(ic)/rhof/omega^2*alpha_f*alpha*beta*H0/Cp*(1-exp(1i*omega*0.03*1d-3))*(-1i)*omega;
end

DC_ubm = real(ubm1(:,1)); DC_vbm = 0*real(vbm1(:,1));
DC_adjBM = 0*real(adjBM(:,1)); DC_adjOW = 0*real(adjOW(1));
DC_ow = real(uow(1)); DC_rw = real(urw(1));

for Null = 1:1  % IFFT process for the original data
    V = zeros(lengths+1,1);
    V(2:nstep+1) = conj(uow(1:nstep));
    V(lengths+1:-1:nstep+2) = uow(1:nstep);
    V(1) = DC_ow;
    ImpulseOW = zeros(lengths+1,2);
    C_ImpulseOW = ImpulseOW;
    ImpulseOW(:,1) = ifft(V*(lengths+1)*df*1d3);
    ImpulseOW(:,2) = ifft(V*(lengths+1)*df*1d3);
    Mow = max(max(abs(real(ImpulseOW))));

    V = zeros(lengths+1,1);
    V(2:nstep+1) = conj(urw(1:nstep));
    V(lengths+1:-1:nstep+2) = urw(1:nstep);
    V(1) = DC_rw;
    ImpulseRW = zeros(lengths+1,2);
    C_ImpulseRW = ImpulseRW;
    ImpulseRW(:,1) = ifft(V*(lengths+1)*df*1d3);
    ImpulseRW(:,2) = ifft(V*(lengths+1)*df*1d3);
    Mrw = max(max(abs(real(ImpulseRW))));

    V = zeros(lengths+1,size(ubm1,1));
    V(2:nstep+1,:) = ubm1(:,1:nstep)';
    V(lengths+1:-1:nstep+2,:) = ubm1(:,1:nstep).';
    V(1,:) = DC_ubm';
    ImpulseBM = ifft(V*(lengths+1)*df*1d3);
    C_ImpulseBM = zeros(lengths+1,size(ubm1,1));
    Mbm = max(max(20*log10(abs(real(ImpulseBM)))));
    
    V = zeros(lengths+1,size(vbm1,1));
    V(2:nstep+1,:) = vbm1(:,1:nstep)';
    V(lengths+1:-1:nstep+2,:) = vbm1(:,1:nstep).';
    V(1,:) = DC_vbm';
    ImpulseVBM = ifft(V*(lengths+1)*df*1d3);

    V = zeros(lengths+1,size(adjBM,1));
    V(2:nstep+1,:) = adjBM(:,1:nstep)';
    V(lengths+1:-1:nstep+2,:) = adjBM(:,1:nstep).';
    V(1,:) = DC_adjBM';
    Impulse_adjBM = ifft(V*(lengths+1)*df*1d3);

    V = zeros(lengths+1,1);
    V(2:nstep+1) = conj(adjOW(1:nstep));
    V(lengths+1:-1:nstep+2,:) = adjOW(1:nstep);
    V(1) = DC_adjOW;
    Impulse_adjOW = ifft(V*(lengths+1)*df*1d3);
end


%% plot

% BM and stapes displacement for acoustic or internal force excitation

figure;
subplot(2,1,2)
n1 = 1;  % 
n2 = 101;  % x position
plot(t,ImpulseOW(:,n1)/max(abs(ImpulseOW(:,n1)))^0,'-c','linewidth',1.); hold on
plot(t,ImpulseBM(:,n2)/max(abs(ImpulseBM(:,n2))),'-k','linewidth',1.); hold on
set(gca,'FontSize',14);
title([num2str((n2-1)*0.025),' mm, max. ',...
    num2str(max(abs(ImpulseBM(:,n2)))*1e1),' mm'],'FontSize',14);
xlabel('msec','FontSize',14);
ylabel('Normalized displacement','FontSize',14);
axis([0 3 -1 1]);



% BM displacement from internal pressure excitation
figure;
n2 = 101;  % x position
plot(t,real(Impulse_adjBM(:,n2))/max(abs(real(Impulse_adjBM(:,n2)))),...
    '-k','linewidth',1); hold off
set(gca,'FontSize',14);
title([num2str((n2-1)*0.025),' mm, max. ',...
    num2str(max(abs(Impulse_adjBM(:,n2)))*1e4),' {\mu}m/s'],'FontSize',14);
xlabel('msec','FontSize',14); ylabel('Normalized velocity','FontSize',14);
axis([-0.1 3 -1 1]);

% The the group delay of each coda
n2 = 161;  % x position
tmpData = real(Impulse_adjBM(:,n2));
t1 = 2.6; % ms
t2 = 3.8; % ms
T1 = round(t1/dt)+1;
T2 = round(t2/dt)+1;
getData = zeros(lengths+1,1);
getData(T1:T2) = tmpData(T1:T2);
figure; plot(t,getData);
fftData = fft(getData);
figure;
subplot(2,1,1)
plot(freq,abs(fftData(1:nsteps)))
subplot(2,1,2)
plot(freq,unwrap(angle(fftData(1:nsteps)))/2/pi)
phData = unwrap(angle(fftData(1:nsteps)))/2/pi;
f1 = 9; % kHz
f2 = 10.4; % kHz
F1 = round(f1/df)+1;
F2 = round(f2/df)+1;
delay = (phData(F2)-phData(F1))/(f2-f1);

% BM velocity from acoustic or internal force excitation
figure;
n2 = 101;  %  x position
plot(t,real(ImpulseVBM(:,n2))/max(abs(real(ImpulseVBM(:,n2)))),...
    '-k','linewidth',1); hold off
set(gca,'FontSize',14);
title([num2str((n2-1)*0.025),' mm, max. ',...
    num2str(max(abs(ImpulseVBM(:,n2)))*1e4),' {\mu}m/s'],'FontSize',14);
xlabel('msec','FontSize',14); ylabel('Normalized velocity','FontSize',14);
axis([-0.1 3 -1 1]);