% University of Rhode Island
% Dynamic Photomechanics Laboratory
 
clear all
close all
 
 
bar_dia = 0.50*25.4/1000;                                                     % the outer diameter of the bar
bar_dia_hollow = 0;                                                           % the inner diameter of the hollow bar
 
bar_e = 200e9;                                                                % the Young's modulus of the bar (MPa)
bar_c = 5010;                                                                 % the elastic wave velocity in the bar (m/s)
inc_bar_a = pi*bar_dia*bar_dia/4;                                             % the cross section area of the incident bar 
tra_bar_a = pi*(bar_dia*bar_dia - bar_dia_hollow*bar_dia_hollow)/4;           % the cross section area of the transmitted bar
 
%__________________________________________________________________________
 
% Input from the oscilloscope
 
disp('Now please input the data filename of channel 1(without extension):');
data_name1='TEK00000';
data_name2='TEK00001';
data_name3='TEK00002';
data_name4='TEK00003';
 
 
disp('Now please input the extension of the files:');
data_extension='csv';
 
eval(['load ',data_name1,'.',data_extension,';'])
eval(['load ',data_name2,'.',data_extension,';'])
eval(['load ',data_name3,'.',data_extension,';'])
eval(['load ',data_name4,'.',data_extension,';'])
disp(' ');
 
%__________________________________________________________________________
 
% Converting the output voltage to microstrains
 
eval(['siganl1= ',data_name1,'(:,2)*1000*1/1.065;'])
eval(['siganl2= ',data_name2,'(:,2)*1000*1/1.065;'])
eval(['siganl3= ',data_name3,'(:,2)*1000*1/1.065;'])
eval(['siganl4= ',data_name4,'(:,2)*1000*1/1.065;'])
 
eval(['dt= ',data_name1,'(2,1)-',data_name1,'(1,1);'])
 
%__________________________________________________________________________
 
% Balacing signal
 
signal1_1 = siganl1(1:500,1);
signal2_1 = siganl2(1:500,1);
signal3_1 = siganl3(1:500,1);
signal4_1 = siganl4(1:500,1);
signal1avg = mean(signal1_1);
signal2avg = mean(signal2_1);
signal3avg = mean(signal3_1);
signal4avg = mean(signal4_1);
siganl1 = siganl1 - signal1avg;
siganl2 = siganl2 - signal2avg;
siganl3 = siganl3 - signal3avg;
siganl4 = siganl4 - signal4avg;
 
%__________________________________________________________________________
 
% Averaging the pulses on the incident and transmission bars
 
inc_pulse_originraw = (siganl1 + siganl2) / 2.0;                            % the signal in incident bar
tra_pulse_originraw = (siganl3 + siganl4) / 2.0;                            % the signal in transmitted bar
 
%__________________________________________________________________________
 
% Filtering the avaraged pulses
 
fn=0.2;
n=2;
 
[b,a] = butter(n, fn );
 
inc_pulse_origin = filtfilt(b,a,inc_pulse_originraw);
 
fn=0.05;
n=2;
 
[bb,a] = butter(n, fn );
 
 
tra_pulse_origin = filtfilt(bb,a,tra_pulse_originraw);   
 
%__________________________________________________________________________
 
figure(1), 
eval(['plot(',data_name1,'(:,1)*1000000,inc_pulse_origin,','''','r','''','),hold on;'])
grid on;
 
figure(2), 
eval(['plot(',data_name1,'(:,1)*1000000,tra_pulse_origin,','''','r','''','),hold on;'])
grid on;
 
pause;
 
eval(['t0= ',data_name1,'(1,1)*1000000;'])
 
beginc=input('please input the time that incident pulse begins:');
n_inc_begin=ceil((beginc-t0)/(dt*1000000))+1;
n_inc_end=ceil((input('please input the time that incident pulse ends:')-t0)/(dt*1000000))+1;
n_ref_begin=ceil((input('please input the time that reflected pulse begins:')-t0)/(dt*1000000))+1;
n_tra_begin=ceil((input('please input the time that transmitted pulse begins:')-t0)/(dt*1000000))+1;
 
n0=n_inc_end-n_inc_begin+1;
 
for i=1:n0;
    t(i,1)=(i-1)*dt;
    inc_pulse(i,1)=inc_pulse_origin(n_inc_begin-1+i);
    ref_pulse(i,1)=inc_pulse_origin(n_ref_begin-1+i);
    tra_pulse(i,1)=tra_pulse_origin(n_tra_begin-1+i);
    
    P_inc(i,1)=((inc_pulse(i,1)+ref_pulse(i,1))/1000000)*bar_e*inc_bar_a;
    P_tra(i,1)=(tra_pulse(i,1)/1000000)*bar_e*tra_bar_a;
end
 
figure(3), 
eval(['plot(',data_name1,'(:,1)*1000000,inc_pulse_origin,','''','r','''','),hold on;'])
eval(['plot(',data_name1,'(:,1)*1000000,tra_pulse_origin,','''','b','''','),hold on;'])
grid on; axis tight;
title('Original pulses','FontName','Timesnewroman','FontSize',12);
xlabel('Time(\mus)','FontName','Timesnewroman','FontSize',12);
ylabel('strain(\mu\epsilon)','FontName','Timesnewroman','FontSize',12);
legend('incident pulse','transmitted pulse','FontName','Timesnewroman','FontSize',12);
 
figure(4) 
plot(t*1000000,inc_pulse,'r'),hold on;
plot(t*1000000,-ref_pulse,'g'),hold on;
plot(t*1000000,tra_pulse,'b'),hold on;
axis tight; grid on;
title('Pulses zoom in','FontName','Timesnewroman','FontSize',12);
xlabel('Time(\mus)','FontName','Timesnewroman','FontSize',12);
ylabel('strain(\mu\epsilon)','FontName','Timesnewroman','FontSize',12);
legend('incident pulse','reflected pulse','transmitted pulse','FontName','Timesnewroman','FontSize',12);
 
figure(5) 
plot(t*1000000,P_inc,'r.'),hold on;
plot(t*1000000,P_tra,'b'),hold on;
title('Force applied on the specimen','FontName','Timesnewroman','FontSize',22);
axis tight; grid on;
xlabel('Time(µs)','FontName','Timesnewroman','FontSize',22);
ylabel('Force(N)','FontName','Timesnewroman','FontSize',22);
legend('Front face','Back face','FontName','Timesnewroman','FontSize',22);