% This programme will generate a graph of superimposed calcualted and experimental velocities.
% Two graphs will be produced, graph in plane 1-2, and graph in plane 1-3.
% 


% Copyright David Paterson, University of Strathclyde, Glasgow, UK. 
% Reuse or reuse with editing is permitted but must cite the author of this work. 
% David A.P. Paterson, "Life Cycle and Ultrasonic Based Non-Destructive Analysis of Recycled and Remanufactured Carbon Fibre Reinforced Plastic Composite",
% Doctoral Thesis, Univeristy of Strathclyde, Glasgow, 2018

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                              %velocities for V-CFRP

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                              
                              
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   sample 1     %%%%%%%%%%%%%%%
Sample1_1_2_VL = [A B C D E F G H etc]; % 1-2 plane for sample 1 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample1_1_2_VT = [I J K L M N O P etc]; % 1-2 plane for sample 1 and transverse velocities from experiment
%Incident Angle = [9,10,11,12,13,14,15,16 etc...]
Sample1_1_3_VL = [A B C D E F G H etc]; % 1-3 plane for sample 3 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample1_1_3_VT = [I J K L M N O P etc]; % 1-3 plane for sample 3 and and transverse velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sample 1    %%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sample 2    %%%%%%%%%%%%%%%
Sample2_1_2_VL = [A B C D E F G H etc]; % 1-2 plane for sample 1 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample2_1_2_VT = [I J K L M N O P etc]; % 1-2 plane for sample 1 and transverse velocities from experiment
%Incident Angle = [9,10,11,12,13,14,15,16 etc...]
Sample2_1_3_VL = [A B C D E F G H etc]; % 1-3 plane for sample 3 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample2_1_3_VT = [I J K L M N O P etc]; % 1-3 plane for sample 3 and and transverse velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sample 2    %%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sample 3    %%%%%%%%%%%%%%%
Sample3_1_2_VL = [A B C D E F G H etc]; % 1-2 plane for sample 1 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample3_1_2_VT = [I J K L M N O P etc]; % 1-2 plane for sample 1 and transverse velocities from experiment
%Incident Angle = [9,10,11,12,13,14,15,16 etc...]
Sample3_1_3_VL = [A B C D E F G H etc]; % 1-3 plane for sample 3 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample3_1_3_VT = [I J K L M N O P etc]; % 1-3 plane for sample 3 and and transverse velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sample 3    %%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sample 4    %%%%%%%%%%%%%%%
Sample4_1_2_VL = [A B C D E F G H etc]; % 1-2 plane for sample 1 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample4_1_2_VT = [I J K L M N O P etc]; % 1-2 plane for sample 1 and transverse velocities from experiment
%Incident Angle = [9,10,11,12,13,14,15,16 etc...]
Sample4_1_3_VL = [A B C D E F G H etc]; % 1-3 plane for sample 3 and longitudinal velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
Sample4_1_3_VT = [I J K L M N O P etc]; % 1-3 plane for sample 3 and and transverse velocities from experiment
%Incident Angle = [0,1,2,3,4,5,6,7,8 etc...]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sample 4    %%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    Average velocities    %%%%%%%%%%%
AvgExp_VL_1_2 = ((Sample1_1_2_VL + Sample2_1_2_VL + Sample3_1_2_VL + Sample4_1_2_VL)/4); %Average velocites from all sample measurements
AvgExp_VT_1_2 = ((Sample1_1_2_VT + Sample2_1_2_VT + Sample3_1_2_VT + Sample4_1_2_VT)/4); %Average velocites from all sample measurements
 
AvgExp_VL_1_3 = ((Sample1_1_3_VL + Sample2_1_3_VL + Sample3_1_3_VL + Sample4_1_3_VL)/4); %Average velocites from all sample measurements
AvgExp_VT_1_3 = ((Sample1_1_3_VT + Sample2_1_3_VT + Sample3_1_3_VT + Sample4_1_3_VT)/4); %Average velocites from all sample measurements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    Average velocities    %%%%%%%%%%%%

% P = ((A/((B))*1000)+(A/((B))*1000)+(A/((B))*1000)+(A/((B))*1000))/4; % density equations (Mass/Volume)) for each samples then divide by total number of samples
% SampleV_in_water = (A+B+C+D)/4;
 
% average constants from average wave velocity
 
% c11 = A;
% c33 = B;
% c13 = C;
% c44 = D;
% c66 = E;
% c12 = c11-(2*(c66));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   Velocity selection    %%%%%%%%%%%%
 
% Exp_1_2_VL  = (Sample1_1_2_VL); % Sample 1 Velocities
% Exp_1_2_VT  = (Sample1_1_2_VT); % Sample 1 Velocities
% Exp_1_3_VL  = (Sample1_1_3_VL); % Sample 1 Velocities
% Exp_1_3_VT  = (Sample1_1_3_VT); % Sample 1 Velocities
% P = (A/((B))*1000);
% SampleV_in_water = (C);
% c = [A B C D E]; %v_CFRP sample 1

%Exp_1_2_VL  = (Sample2_1_2_VL); % Sample 2 Velocities
%Exp_1_2_VT  = (Sample2_1_2_VT); % Sample 2 Velocities
%Exp_1_3_VL  = (Sample2_1_3_VL); % Sample 2 Velocities
%Exp_1_3_VT  = (Sample2_1_3_VT); % Sample 2 Velocities
% P = (A/((B))*1000);
% SampleV_in_water = (C);
%c = [A B C D E]; %v_CFRP sample 2

% Exp_1_2_VL  = (Sample3_1_2_VL); % Sample 3 Velocities
% Exp_1_2_VT  = (Sample3_1_2_VT); % Sample 3 Velocities
% Exp_1_3_VL  = (Sample3_1_3_VL); % Sample 3 Velocities
% Exp_1_3_VT  = (Sample3_1_3_VT); % Sample 3 Velocities
% P = (A/((B))*1000);
% SampleV_in_water = (C);
%c = [A B C D E]; %v_CFRP sample 3

% Exp_1_2_VL  = (Sample4_1_2_VL); % Sample 4 Velocities
% Exp_1_2_VT  = (Sample4_1_2_VT); % Sample 4 Velocities
% Exp_1_3_VL  = (Sample4_1_3_VL); % Sample 4 Velocities
% Exp_1_3_VT  = (Sample4_1_3_VT); % Sample 4 Velocities
% P = (A/((B))*1000);
% SampleV_in_water = (C);
%c = [A B C D E]; %v_CFRP sample 4

% c11 = c(1); % assigns c(1) to c11
% c33 = c(2); % assigns c(1) to c33
% c44 = c(3); % assigns c(1) to c44
% c13 = c(4); % assigns c(1) to c13
% c66 = c(5); % assigns c(1) to c66
% c12 = c11-(2*c66); % calculates c12

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    Velocity selection    %%%%%%%%%%%%

% These graphs need to match the array sizes of the velocities recorded above in order for programme to execute correctly 

NewIncident = 0:45; % sets up an array of variables needed for graphing purposes. This has 46 entries

% =========  assigns the experimental refraction velocity to an array incident angle  ============================

IncExpVqt = zeros(1,46); % sets up the average velocites into an array for graphing purposes - first array entry is 0 degrees incident
IncExpVqt(17)= Exp_1_2_VT(1); % this gives Exp_1_2_VT  at incident angle of 16 degrees
IncExpVqt(18)= Exp_1_2_VT(2); % this gives Exp_1_2_VT  at incident angle of 17 degrees
IncExpVqt(19)= Exp_1_2_VT(3); % this gives Exp_1_2_VT  at incident angle of 18 degrees
IncExpVqt(20)= Exp_1_2_VT(4); % this gives Exp_1_2_VT  at incident angle of 19 degrees
IncExpVqt(21)= Exp_1_2_VT(5); % this gives Exp_1_2_VT  at incident angle of 20 degrees
IncExpVqt(22)= Exp_1_2_VT(6); % this gives Exp_1_2_VT  at incident angle of 21 degrees
IncExpVqt(25)= Exp_1_2_VT(7); % this gives Exp_1_2_VT  at incident angle of 24 degrees
IncExpVqt(28)= Exp_1_2_VT(8); % this gives Exp_1_2_VT  at incident angle of 27 degrees
IncExpVqt(31)= Exp_1_2_VT(9); % this gives Exp_1_2_VT  at incident angle of 30 degrees
IncExpVqt(34)= Exp_1_2_VT(10); % this gives Exp_1_2_VT  at incident angle of 33 degrees
IncExpVqt(37)= Exp_1_2_VT(11); % this gives Exp_1_2_VT  at incident angle of 36 degrees

IncExpVl = zeros(1,46); % sets up the average velocites into an array for graphing purposes - first array entry is 0 degrees incident
IncExpVl(1)= Exp_1_2_VL(1); % this gives Exp_1_2_VL  at incident angle of 0 degrees
IncExpVl(4)= Exp_1_2_VL(2); % this gives Exp_1_2_VL  at incident angle of 3 degrees
IncExpVl(7)= Exp_1_2_VL(3); % this gives Exp_1_2_VL  at incident angle of 6 degrees
IncExpVl(10)= Exp_1_2_VL(4); % this gives Exp_1_2_VL  at incident angle of 9 degrees
IncExpVl(13)= Exp_1_2_VL(5); % this gives Exp_1_2_VL  at incident angle of 12 degrees
IncExpVl(15)= Exp_1_2_VL(6); % this gives Exp_1_2_VL  at incident angle of 14 degrees
IncExpVl(16)= Exp_1_2_VL(7); % this gives Exp_1_2_VL  at incident angle of 15 degrees

Calc_VT = sqrt(c66/P); % calculated velocity using elastic constants value
CalcVqt = zeros(1,46);

CalcVqt(17)= Calc_VT; % this gives Calc_VT at incident angle of 16 degrees
CalcVqt(18)= Calc_VT; % this gives Calc_VT at incident angle of 17 degrees
CalcVqt(19)= Calc_VT; % this gives Calc_VT at incident angle of 18 degrees
CalcVqt(20)= Calc_VT; % this gives Calc_VT at incident angle of 19 degrees
CalcVqt(21)= Calc_VT; % this gives Calc_VT at incident angle of 20 degrees
CalcVqt(22)= Calc_VT; % this gives Calc_VT at incident angle of 21 degrees
CalcVqt(25)= Calc_VT; % this gives Calc_VT at incident angle of 24 degrees
CalcVqt(28)= Calc_VT; % this gives Calc_VT at incident angle of 27 degrees
CalcVqt(31)= Calc_VT; % this gives Calc_VT at incident angle of 30 degrees
CalcVqt(34)= Calc_VT; % this gives Calc_VT at incident angle of 33 degrees
CalcVqt(37)= Calc_VT; % this gives Calc_VT at incident angle of 36 degrees

Calc_VL = sqrt(c11/P); % calculated velocity using elastic constants value

CalcVl = zeros(1,46); % sets up the average velocites into an array for graphing purposes - first array entry is 0 degrees incident
CalcVl(1)= Calc_VL; % this gives Calc_VL at incident angle of 0 degrees
CalcVl(4)= Calc_VL; % this gives Calc_VL at incident angle of 3 degrees
CalcVl(7)= Calc_VL; % this gives Calc_VL at incident angle of 6 degrees
CalcVl(10)= Calc_VL; % this gives Calc_VL at incident angle of 9 degrees
CalcVl(13)= Calc_VL; % this gives Calc_VL at incident angle of 12 degrees
CalcVl(15)= Calc_VL; % this gives Calc_VL at incident angle of 14 degrees
CalcVl(16)= Calc_VL; % this gives Calc_VL at incident angle of 15 degrees

%========================  velocity Graph  ===================================

figure(1)
scatter(NewIncident,IncExpVqt,'filled') % data from experimental 
hold on
scatter(NewIncident,IncExpVl,'filled') % data from experimental
scatter(NewIncident,CalcVqt,'filled') % data from experimental
scatter(NewIncident,CalcVl,'filled') % data from experimental
axis([0 45 1 3500]) % Sets the axis values for both x and y
grid on
xlabel('Incident Angle (deg)'); % labels the x axis
ylabel('Wave velocity (m/s)'); % Labels the y axis
title('Experimental\Calculated velocity against incident angle for plane 1-2 v-CFRP sample average') % Labels the graph

%========================  velocity Graph  ===================================

% These graphs need to match the array sizes of the velocities recorded above in order for programme to execute correctly 

% =========  assigns the experimental refraction velocity to an array incident angle  ============================

IncExpVqt = zeros(1,46); % sets up the average velocites into an array for graphing purposes - first array entry is 0 degrees incident
IncExpVqt(6)= Exp_1_3_VT (1); % this gives Exp_1_3_VT  at incident angle of 5 degrees
IncExpVqt(7)= Exp_1_3_VT (2); % this gives Exp_1_3_VT  at incident angle of 6 degrees
IncExpVqt(9)= Exp_1_3_VT (3); % this gives Exp_1_3_VT  at incident angle of 8 degrees
IncExpVqt(11)= Exp_1_3_VT (4); % this gives Exp_1_3_VT  at incident angle of 10 degrees
IncExpVqt(13)= Exp_1_3_VT (5); % this gives Exp_1_3_VT  at incident angle of 12 degrees
IncExpVqt(15)= Exp_1_3_VT (6); % this gives Exp_1_3_VT  at incident angle of 14 degrees
IncExpVqt(17)= Exp_1_3_VT (7); % this gives Exp_1_3_VT  at incident angle of 16 degrees
IncExpVqt(19)= Exp_1_3_VT (8); % this gives Exp_1_3_VT  at incident angle of 18 degrees
IncExpVqt(21)= Exp_1_3_VT (9); % this gives Exp_1_3_VT  at incident angle of 20 degrees
IncExpVqt(23)= Exp_1_3_VT (10); % this gives Exp_1_3_VT  at incident angle of 22 degrees
IncExpVqt(25)= Exp_1_3_VT (11); % this gives Exp_1_3_VT  at incident angle of 24 degrees
IncExpVqt(27)= Exp_1_3_VT (12); % this gives Exp_1_3_VT  at incident angle of 26 degrees
IncExpVqt(29)= Exp_1_3_VT (13); % this gives Exp_1_3_VT  at incident angle of 28 degrees
IncExpVqt(31)= Exp_1_3_VT (14); % this gives Exp_1_3_VT  at incident angle of 30 degrees
IncExpVqt(33)= Exp_1_3_VT (15); % this gives Exp_1_3_VT  at incident angle of 32 degrees

IncExpVl = zeros(1,46); % sets up the average velocites into an array for graphing purposes - first array entry is 0 degrees incident
IncExpVl(1)= Exp_1_3_VL (1); % this gives Exp_1_3_VT L at incident angle of 0 degrees
IncExpVl(3)= Exp_1_3_VL (2); % this gives Exp_1_3_VT L at incident angle of 2 degrees
IncExpVl(4)= Exp_1_3_VL (3); % this gives Exp_1_3_VT L at incident angle of 3 degrees
IncExpVl(5)= Exp_1_3_VL (4); % this gives Exp_1_3_VT L at incident angle of 4 degrees

%Outputs the graph of wave velocity

%=================================== Refraction angles transverse ===========================================================================

while counter == 1 
    
    Incident = [A B C D E F G H I J K etc]; % sets up the Incident angles
    sindrefract(counter2) = ((sind(Incident(counter2)))*Exp_VT_1_3(counter2))/(SampleV_in_water); % calculates the sinde of refraction angle
    Refract(counter2) = asind(sindrefract(counter2)); % calcualtes the refraction anlge (this is the anlge the way travels through sample in)
   
    counter2 = counter2+1; % increases counter2 variable to allow for  moves to create new refraction angle.
    
    if counter2 == A % (limited here by number of incident angles)
        counter = 0; % sets the while loop counter value and exists while loop
    else
        counter = 1; % sets the while loop counter value
    end % end of if else statement
   
end % end of while loop 

%=================================== Refraction angles transverse ===========================================================================

%====================================  Refraction angles longitudinal  ======================================================================

counter = 1; % resets the counter
counter2 = 1; % resets the counter

while counter == 1 
    
    IncidentL = [A B C D E F G H I J K etc]; % sets up the Incident angles
    sindrefractL(counter2) = ((sind(IncidentL(counter2)))*Exp_VL_1_3(counter2))/(SampleV_in_water); % calculates the sinde of refraction angle
    RefractL(counter2) = asind(sindrefractL(counter2)); % calcualtes the refraction anlge (this is the anlge the way travels through sample in)
  
    counter2 = counter2+1; % increases counter2 variable to allow for  moves to create new refraction angle.
    
    if counter2 == A %(limited here by number of incident angles)
        counter = 0; % sets the while loop counter value and exists while loop
    else
        counter = 1; % sets the while loop counter value
    end % end of if else statement
    
end % end of while loop 

% ====================================  Refraction angles longitudinal  ======================================================================
 
counter = 1; % resets the variable to allow it to be used later on
counter2 = 1; % resets the variable to allow it to be used later on

% ===================== Calculated Transverse Velocities  ==============================
while counter == 1 

    A(counter2) = ((c11*(cosd(Refract(counter2))^2))+(c33*(sind(Refract(counter2))^2))+c44);
    B(counter2) = ((c11*(cosd(Refract(counter2))^2))+(c44*(sind(Refract(counter2))^2)))*((c44*(cosd(Refract(counter2))^2))+(c33*(sind(Refract(counter2))^2)))-(((c13+c44)^2)*(sind(Refract(counter2))^2)*(cosd(Refract(counter2))^2));
    Vqt(counter2) = sqrt((A(counter2)-sqrt((A(counter2)^2)-(4*B(counter2))))/(2*P));

    counter2 = counter2+1; % increases counter2 variable to allow for  moves to create new refraction angle.
    
    if counter2 == A % (limited here by number of refract angles)
        counter = 0; % sets the while loop counter value and exists while loop
    else
        counter = 1; % sets the while loop counter value
    end % end of if else statement
    
    %S = sum(Vqt);
end % end of while loop 
   
% ====================================  Calculated Transverse Velocities  ======================================================================

counter = 1;
counter2 = 1;

% ====================================  calculated Longitudinal Velocities ======================================================================

while counter == 1 
    
    C(counter2) = ((c11*(cosd(RefractL(counter2))^2))+(c33*(sind(RefractL(counter2))^2))+c44);
    D(counter2) = ((c11*(cosd(RefractL(counter2))^2))+(c44*(sind(RefractL(counter2))^2)))*((c44*(cosd(RefractL(counter2))^2))+(c33*(sind(RefractL(counter2))^2)))-(((c13+c44)^2)*(sind(RefractL(counter2))^2)*(cosd(RefractL(counter2))^2));
    Vl(counter2) = sqrt((C(counter2) + sqrt((C(counter2)^2)-(4*D(counter2))))/(2*P));
   
    counter2 = counter2+1; % increases counter2 variable to allow for  moves to create new refraction angle.
    
    if counter2 == A % (limited here by number of refractL angles)
        counter = 0; % sets the while loop counter value and exists while loop
    else
        counter = 1; % sets the while loop counter value
    end % end of if else statement
  
end % end of while loop 

% ====================================  Calculated Longitudinal Velocities ======================================================================

% =========  assigns the experimental refraction velocity to an array incident angle  ============================

% These graphs need to match the array sizes of the velocities recorded above in order for programme to execute correctly 

IncCalcVqt = zeros(1,46); % sets up the average velocites into an array for graphing purposes - first array entry is 0 degrees incident
IncCalcVqt(6) = Vqt(1); % this gives Vqt at incident angle of 5 degrees
IncCalcVqt(7) = Vqt(2); % this gives Vqt at incident angle of 6 degrees
IncCalcVqt(9) = Vqt(3); % this gives Vqt at incident angle of 8 degrees
IncCalcVqt(11) = Vqt(4); % this gives Vqt at incident angle of 10 degrees
IncCalcVqt(13) = Vqt(5); % this gives Vqt at incident angle of 12 degrees
IncCalcVqt(15) = Vqt(6); % this gives Vqt at incident angle of 14 degrees
IncCalcVqt(17) = Vqt(7); % this gives Vqt at incident angle of 16 degrees
IncCalcVqt(19) = Vqt(8); % this gives Vqt at incident angle of 18 degrees
IncCalcVqt(21) = Vqt(9); % this gives Vqt at incident angle of 20 degrees
IncCalcVqt(23) = Vqt(10); % this gives Vqt at incident angle of 22 degrees
IncCalcVqt(25) = Vqt(11); % this gives Vqt at incident angle of 24 degrees
IncCalcVqt(27) = Vqt(12); % this gives Vqt at incident angle of 26 degrees
IncCalcVqt(29) = Vqt(13); % this gives Vqt at incident angle of 28 degrees
IncCalcVqt(31) = Vqt(14); % this gives Vqt at incident angle of 30 degrees
IncCalcVqt(33) = Vqt(15); % this gives Vqt at incident angle of 32 degrees

IncCalcVl = zeros(1,46); % sets up the average velocites into an array for graphing purposes - first array entry is 0 degrees incident
IncCalcVl(1) = Vl(1); % this gives Vl at incident angle of 0 degrees
IncCalcVl(3) = Vl(2); % this gives Vl at incident angle of 2 degrees 
IncCalcVl(4) = Vl(3); % this gives Vl at incident angle of 3 degrees
IncCalcVl(5) = Vl(4); % this gives Vl at incident angle of 4 degrees

%======================== velocity Graph  ===================================

figure(2)
scatter(NewIncident,IncExpVqt,'filled') % data from experimental 
hold on
scatter(NewIncident,IncExpVl,'filled') % data from experimental
scatter(NewIncident,IncCalcVqt,'filled') % data from experimental
scatter(NewIncident,IncCalcVl,'filled') % data from experimental
axis([0 45 1 3500]) % Sets the axis values for both x and y
grid on
xlabel('Incident Angle (deg)'); % labels the x axis
ylabel('Wave velocity (m/s)'); % Labels the y axis
title('Experimental\Calculated velocity against incident angle for plane 1-3 v-CFRP sample average') % Labels the graph

%======================== velocity Graph  ===================================

% DIFFERENCE CALCULATIONS. This section of code calcuates % difference between the difference
% vetw

counter = 1;
counter2 = 1;

while counter == 1
PD_VT_1_2(counter2) = (abs(Calc_VT-Exp_1_2_VT(counter2)))/((Calc_VT+Exp_1_2_VT(counter2))/2)*100;

counter2 = counter2+1;

    if counter2 == A % (limited here by number of wave measurments)
        counter = 0;
    else
        counter = 1;
    end
end

counter = 1;
counter2 = 1;

while counter == 1    
PD_VL_1_2(counter2) = (abs(Calc_VL-Exp_1_2_VL(counter2)))/((Calc_VL+Exp_1_2_VL(counter2))/2)*100;

counter2 = counter2+1;

    if counter2 == A % (limited here by number of wave measurments)
        counter = 0;
    else
        counter = 1;
    end
end

counter = 1;
counter2 = 1;

while counter == 1
PD_VT_1_3(counter2) = (abs(Vqt(counter2)-Exp_1_3_VT(counter2)))/((Vqt(counter2)+Exp_1_3_VT(counter2))/2)*100;

counter2 = counter2+1;

    if counter2 == A % (limited here by number of wave measurments)
        counter = 0;
    else
        counter = 1;
    end
end

counter = 1;
counter2 = 1;

while counter == 1    
PD_VL_1_3(counter2) = (abs(Vl(counter2)-Exp_1_3_VL(counter2)))/((Vl(counter2)+Exp_1_3_VL(counter2))/2)*100;

counter2 = counter2+1;

    if counter2 == A % (limited here by number of wave measurments)
        counter = 0;
    else
        counter = 1;
    end
end

Exp_1_3_VT % example of how to select % difference figures
Vqt 
PD_VT_1_3