% University of Rhode Island
% Dynamic Photomechanics Laboratory
% This code was developed for SHPB analysis

clear all
close all
 
% Bar and specimen parameters
 
bar_dia = 12.7/1000;                                                         % the outer diameter of the bar
bar_dia_hollow = 0/1000;                                                     % the inner diameter of the hollow bar
 
spe_dia=input ('enter specimen diameter in inches: ');                       % the diameter of the specimen
spe_l=input ('enter specimen thickness in inches: ');                        % the thickness of the specimen
 
spe_dia = spe_dia*25.4/1000;                                                 % change the diameter of the specimen into m
spe_l = spe_l*25.4/1000;                                                     % change the thickness of the specimen into m
 
spe_ai = pi*spe_dia*spe_dia/4;                                               % the initial cross section area of the specimen (m^2) 
bar_e = 200e3;                                                               % 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
con = bar_c/spe_l;                                                           % constant in the eqations 
%__________________________________________________________________________
 
% 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(['time_origin= ',data_name1,'(:,1);']);
 
dt=time_origin(2,1)-time_origin(1,1)                                        % time between two points of the data (s)
 
%__________________________________________________________________________
 
% 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_origin = (siganl1 + siganl2) / 2.0;                               % the signal in incident bar
tra_pulse_origin = (siganl3 + siganl4) / 2.0;                               % the signal in transmitted bar
 
%__________________________________________________________________________
 
% Filtering the avaraged pulses
 
fn=0.2;
n=2;
 
[bb,a] = butter(n, fn );
 
inc_pulse = filtfilt(bb,a,inc_pulse_origin);
 
fn=0.05;
n=2;
 
[bbc,a] = butter(n, fn );
 
tra_pulse = filtfilt(bbc,a,tra_pulse_origin);   
 
%__________________________________________________________________________
 
 
figure(1), 
plot(time_origin*1e6,inc_pulse_origin,'r')
title('please find the time that the incident and reflected pulse begin and ends')
grid on;
 
figure(2), 
plot(time_origin*1e6,tra_pulse_origin,'r')
title('please find the time that the transmitted pulse begin and ends')
grid on;
 
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;
 
pulse_length = n_inc_end - n_inc_begin + 1;                                 % Pulse length
 
for i = 1: pulse_length
    k1 = n_ref_begin + i -1;
    k2 = n_tra_begin + i -1;
    pulse_time(i,1) = dt*(i-1);
    refl(i,1) = inc_pulse(k1);
    ref_data(i,1)=dt*(i-1);
    ref_data(i,2)=inc_pulse(k1);
    
    trans(i,1) = tra_pulse(k2);
    tra_data(i,1)=dt*(i-1);
    tra_data(i,2)=tra_pulse(k2);
end
 
%__________________________________________________________________________
 
 
rarea(1)=0;
Rfarea(1)=0;
for n=2:pulse_length,
    rarea(n)=(refl(n-1)+refl(n))*(0.5*dt);
    Rfarea(n)=Rfarea(n-1)+rarea(n);
end
 
% Integral value of the transmitted pulse
tarea(1)=0;
TRarea(1)=0;
for n=2:pulse_length,
    tarea(n)=(trans(n-1)+trans(n))*(0.5*dt);
    TRarea(n)=TRarea(n-1)+tarea(n);
end
 
% Calsulate the true strain & strain rate
for nn=1:pulse_length
    estrain(nn,1)=-con*(((tra_bar_a/inc_bar_a)-1)*TRarea(nn)-2*Rfarea(nn))/1e6;  % engineering strain
    srate(nn,1)=-con*(((tra_bar_a/inc_bar_a)-1)*trans(nn)-2*refl(nn))/1e6;       % engineering rate
    estress(nn,1) = -bar_e*(tra_bar_a/spe_ai)*trans(nn)/1e6;                     % engineering stress
    tstrain(nn,1) = -log(1-estrain(nn,1));                                       % true strain
    tstress(nn,1) = estress(nn,1)*(1-estrain(nn,1));                             % true stress
    true_stress_strain(nn,1)=  tstrain(nn,1);
    true_stress_strain(nn,2)= tstress(nn,1);
    e_stress_strain(nn,1)=  estrain(nn,1);
    e_stress_strain(nn,2)= estress(nn,1);
end
 
%__________________________________________________________________________
 
% PLOT THE CURVES IN MATLAB FOR A QUICK PERUSAL!
 
figure(3)
plot(estrain*100,estress)
xlabel('Engineering Strain (%)','FontName','Timesnewroman','FontSize',12)
ylabel('Engineering Stress (MPa)','FontName','Timesnewroman','FontSize',12)
grid on;
 
figure(4)
plot(tstrain*100,tstress)
xlabel('True Strain (%)','FontName','Timesnewroman','FontSize',12)
ylabel('True Stress (MPa)','FontName','Timesnewroman','FontSize',12)
grid on;
title('Pick first point of two points to calculate the slope, right button to continue');
[x1,y1] = ginput(1);
title('Pick second point of two points to calculate the slope, right button to continue');
[x2,y2] = ginput(1);
slope = (y2 - y1)*100/(x2 - x1)
 
title('Press right mouse button to continue, any other to redo')
    [junkx,junky,click]=ginput(1);
 
figure(5)
plot(pulse_time*1e6,tstrain*100);
xlabel('Time(\mus)','FontName','Timesnewroman','FontSize',22)
ylabel('True Strain (%)','FontName','Timesnewroman','FontSize',22)
grid on;
%title('Pick first point of two points to calculate the strain rate, right button to quit');
[x1,y1] = ginput(1);
%title('Pick second point of two points to calculate the strain rate, right button to quit');
[x2,y2] = ginput(1);
strainrate = (y2 - y1)/((x2 - x1)*100)*1e6
 
figure(6)
plot(pulse_time*1e6,srate);
xlabel('Time(\mus)','FontName','Timesnewroman','FontSize',22)
ylabel('Strain Rate (s^-1)','FontName','Timesnewroman','FontSize',22)
grid on;
 
%__________________________________________________________________________
 
% Saving the data
 
save true.txt true_stress_strain -ascii
save eng.txt e_stress_strain -ascii