clc
clear all
close all

%% Data loading

load demo1.mat

data(isnan(data)) = 0; %Remove the NAN values

t = data(:,1);
p = data(:,2);
v = data(:,3);

rho = 1e3; %kg/m3
%% Measures of an ultrasound burst (10p)

%a)
figure

subplot(211)
plot(t*1e6,v)
title('Driving signal')
xlabel('time, us')
ylabel('voltage,V')

subplot(212)
plot(t*1e6,p)
title('Pressure')
xlabel('time, us')
ylabel('voltage,V')

%b)
PRP = (1050 - 50)*1e-6; % s
%c)
PRF = 1/PRP ; % Hz

%d)
dt = t(2)-t(1);
fs = 1/dt;
L = length(p);

P = abs(fft(p-mean(p))); 
P = P(1:L/2+1);
f = [1:L]*fs/L;
f = f(1:L/2+1);
P  = 20*log10(P/max(P));

figure
plot(f,P) 
title('Amplitude Spectrum of p(t))')
xlabel('f (Hz)')
ylabel('|P(f)|db')

% e)
b_6dB = f(find(P>-6));
b_6dB = [b_6dB(1), b_6dB(end)];

%f)
peak_f = f(find(P==max(P)));    

%g)
center_f = sum(b_6dB)/2;

%h) max_freq ~ 3MHz;

%% Speed of sound in water (5p)
close all

figure
plot(t*1e6,p,t*1e6,v)
xlabel('time ,us')
ylabel('vol tage, V')
legend('hydrophone signal','driving signal' )

%a)
t_travel = 50.26*1e-6;

%b)
R = 75*1e-3 ; %m
c = R/t_travel; %m/s

%c) from 'Measuring the speed of sound in distilled water, S. S. Sekoyan': https://link.springer.com/article/10.1007\%2FBF00991734

T1 = 20.619;  %celsius
T2 = 23.935 ;

v1 = 1484.65 ; %m/s
v2 = 1493.19 ; %m/s

% y = m*x+q
m = (T2-T1)/(v2-v1);
q = T1-m*v1;
T_w = m*c + q;

%% Pressure and intensity measures of an ultrasound burst (15p)
close all
%a)
s = 46e-3; % sensitivity, [V/MPa]
p2 =p/s*1e6; %Pa

figure
plot(t*1e6,p2*1e-6)
title('Pressure')
xlabel('time, us')
ylabel('pressure, MPa')


%b)
PPP = max(p2);

%c)
PNP = min(p2);

%d)
pulse_start_indx  = 32450; %from the graph
pulse_end_indx  = 32650;

I_int = cumtrapz(t(pulse_start_indx:pulse_end_indx),p2(pulse_start_indx:pulse_end_indx).^2); %intensity integral over the pulse duration

figure
plot(t(pulse_start_indx:pulse_end_indx)*1e6-t (pulse_start_indx)*1e6, I_int/max(I_int))
xlabel('time, us')
ylabel('Normalized integrated intensity')

t1 = find(I_int>0.1*I_int(end),1);
t2 = find(I_int>0.9*I_int(end),1);

PD = (t2-t1)*dt*1.25;

%e)
I_PPP =  (PPP)^2/( rho *c )/1e4;  %W/cm2

%f)
I_PNP =  (PNP)^2/( rho *c )/1e4;  %W/cm2

%g) 
I_sppa = I_int(end)/( rho *c )/PD/1e4;  %W/cm2

%h) 
I_spta = I_int(end)/( rho *c )/PRP/1e4;  %W/cm2

%i)
DC_driv_sig = (31849-31795)/(44816-31795)*100;

%j)
DC_pressure = PD/PRP*100;