%Elastic Constant Graphs
% This programme will call the Least_squares_programme and produce all the
% graps from the elastic constant data calculated using the least squares
% minimization technique. This is an iterative programme and so the graphs
% produced will be using various elastic constants. 50 graphs produced in
% total, and the minized value of least square tecnhique will detail the
% best result - the most accurat constant values.
% 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
v_Least_Squares; % the script file with elastic constant and experimental velocity data will be loaded. v_Least_Squares must be in method 1.
% ======================== Counters ==========================================================================
array_counter = 1; % counter needed to allow velocity to be created from elastic constants
counter = 1; % counter variable required for later on in programme
counter2 = 1; % counter variable required for later on in programme
counter3 = 1; % counter variable required for graphing purposes
NewIncident = 0:45; % sets up an array of variables needed for graphing purposes
constcounter = 1;
% ======================== Counters ==========================================================================
% ============================= Main body of programme =======================================================
while constcounter <51
%============ sets up the elastic constants values from the array data
c11 = c_min_array(array_counter,6); % creates c11 from array data
c33 = c_min_array(array_counter,7); % creates c33 from array data
c44 = c_min_array(array_counter,8); % creates c44 from array data
c13 = c_min_array(array_counter,9); % creates c13 from array data
c66 = c_min_array(array_counter,10); % creates c66 from array data
%============ sets up the elastic constants values from the array data
% ===================== 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 ======================================================================
% These graphs need to match the array sizes of the velocities recorded above (from v_Least_Squares) 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(7)= Exp_VT_1_3(1); % this gives Exp_VT at incident angle of 6 degrees
IncExpVqt(9)= Exp_VT_1_3(2); % this gives Exp_VT at incident angle of 8 degrees
IncExpVqt(11)= Exp_VT_1_3(3); % this gives Exp_VT at incident angle of 10 degrees
IncExpVqt(13)= Exp_VT_1_3(4); % this gives Exp_VT at incident angle of 12 degrees
IncExpVqt(15)= Exp_VT_1_3(5); % this gives Exp_VT at incident angle of 14 degrees
IncExpVqt(17)= Exp_VT_1_3(6); % this gives Exp_VT at incident angle of 16 degrees
IncExpVqt(19)= Exp_VT_1_3(7); % this gives Exp_VT at incident angle of 18 degrees
IncExpVqt(21)= Exp_VT_1_3(8); % this gives Exp_VT at incident angle of 20 degrees
IncExpVqt(23)= Exp_VT_1_3(9); % this gives Exp_VT at incident angle of 22 degrees
IncExpVqt(25)= Exp_VT_1_3(10); % this gives Exp_VT at incident angle of 24 degrees
IncExpVqt(27)= Exp_VT_1_3(11); % this gives Exp_VT at incident angle of 26 degrees
IncExpVqt(29)= Exp_VT_1_3(12); % this gives Exp_VT at incident angle of 28 degrees
IncExpVqt(31)= Exp_VT_1_3(13); % this gives Exp_VT at incident angle of 30 degrees
IncExpVqt(33)= Exp_VT_1_3(14); % this gives Exp_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_VL_1_3(1); % this gives Exp_VTL at incident angle of 0 degrees
IncExpVl(3)= Exp_VL_1_3(2); % this gives Exp_VTL at incident angle of 2 degrees
IncExpVl(4)= Exp_VL_1_3(3); % this gives Exp_VTL at incident angle of 3 degrees
IncExpVl(5)= Exp_VL_1_3(4); % this gives Exp_VTL at incident angle of 4 degrees
% ========= assigns the refaction velocity at a array to corrospoding incident angle ============================
% ========= assigns the experimental refraction velocity to an array incident angle ============================
% These graphs need to match the array sizes of the velocities recorded above (from v_Least_Squares) 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(7) = Vqt(1); % this gives Vqt at incident angle of 6 degrees
IncCalcVqt(9) = Vqt(2); % this gives Vqt at incident angle of 8 degrees
IncCalcVqt(11) = Vqt(3); % this gives Vqt at incident angle of 10 degrees
IncCalcVqt(13) = Vqt(4); % this gives Vqt at incident angle of 12 degrees
IncCalcVqt(15) = Vqt(5); % this gives Vqt at incident angle of 14 degrees
IncCalcVqt(17) = Vqt(6); % this gives Vqt at incident angle of 16 degrees
IncCalcVqt(19) = Vqt(7); % this gives Vqt at incident angle of 18 degrees
IncCalcVqt(21) = Vqt(8); % this gives Vqt at incident angle of 20 degrees
IncCalcVqt(23) = Vqt(9); % this gives Vqt at incident angle of 22 degrees
IncCalcVqt(25) = Vqt(10); % this gives Vqt at incident angle of 24 degrees
IncCalcVqt(27) = Vqt(11); % this gives Vqt at incident angle of 26 degrees
IncCalcVqt(29) = Vqt(12); % this gives Vqt at incident angle of 28 degrees
IncCalcVqt(31) = Vqt(13); % this gives Vqt at incident angle of 30 degrees
IncCalcVqt(33) = Vqt(14); % 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
% ========= assigns the experimental refraction velocity to an array incident angle ============================
%======================== Exp and Calc velocity Graph ===================================
%figure (counter3) % graph of calculated transverse velocity
figure;
scatter(NewIncident,IncExpVqt,'filled') % data from experimental
hold on
scatter(NewIncident,IncExpVl,'filled') % data from experimental
axis([0 45 100 4600]) % Sets the axis values for both x and y
grid on
scatter(NewIncident,IncCalcVqt,'filled') % data from calculated
scatter(NewIncident,IncCalcVl,'filled') % data from calculated
xlabel('Incident Angle (deg)'); % labels the x axis
ylabel('Wave velocity (m/s)'); % Labels the y axis
title('Experimental/Calculated velocities against incident angle for plane 1-3 v-CFRP average') % Labels the graph
%======================== Exp and Calc velocity Graph ===================================
array_counter = array_counter + 1; % increases array_counter so next elastic constant values can be loaded
counter3 = counter3 + 1; % increases array_counter so next elastic constant values can be loaded
constcounter = constcounter + 1; % resets the counter to allow the velocities to be calculated again
counter = 1; % resets the counter to allow the velocities to be calculated again
counter2 = 1; % resets the counter to allow the velocities to be calculated again
end % end of while loop
c_min_array % outputs the array to user, in this case the user is able to see which
%======================== Additinal user Graphs =========================================================
%%%%%%%%%%%%%%%%%%%%%% Seperated velocities %%%%%%%%%%%%%%%%%%%%%%%%%
% figure(1) % graph of experimental transverse velocity
% scatter(NewIncident,IncExpVqt)
% grid on
% xlabel('Incident Angle (deg)'); % labels the x axis
% ylabel('Wave velocity (m/s)'); % Labels the y axis
% title('Experimental Transverse Wave velocity against incident angle for plane 1-3 virgin CFRP') % Labels the graph
% axis([0 45 1000 2600])
% close (figure(1))
%
% figure(2) % graph of experimental longitudinal velocity
% scatter(NewIncident,IncExpVl)
% grid on
% xlabel('Incident Angle (deg)'); % labels the x axis
% ylabel('Wave Amplitude (V pk/pk)'); % Labels the y axis
% title('ExperimentalLongitudinal Wave velocity against incident angle for plane 1-3 virgin CFRP') % Labels the graph
% axis([0 45 1000 2600])
% close (figure(2))
%%%%%%%%%%%%%%%%%%%%%% Seperated velocities %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% long and Transverse Calc and Exp velocities %%%%%%%%%%%%
% figure (3) % graph of calculated transverse velocity
% scatter(NewIncident,IncExpVqt,'filled')
% hold on
% scatter(NewIncident,IncExpVl,'filled')
% axis([0 45 100 3600]) % Sets the axis values for both x and y
% grid on
% xlabel('Incident Angle (deg)'); % labels the x axis
% ylabel('Wave velocity'); % Labels the y axis
% title('Experimental Wave velocities against incident angle for plane 1-3 virgin CFRP') % Labels the graph
%
% figure (4) % graph of calculated longtudinal velocity
% scatter(NewIncident,IncCalcVqt,'filled')
% hold on
% scatter(NewIncident,IncCalcVl,'filled')
% axis([0 45 100 3600]) % Sets the axis values for both x and y
% grid on
% xlabel('Incident Angle (deg)'); % labels the x axis
% ylabel('Wave velocity'); % Labels the y axis
% title('Calculated Wave velocities against incident angle for plane 1-3 virgin CFRP') % Labels the graph
%%%%%%%%%%%% long and Transverse Calc and Exp velocities %%%%%%%%%%%%