close all
clear all
clc
 
tic
for k=1:35
 
    eval(sprintf('load C:\\School_Work\\ESM_4714\\Hoop_Strain_Mode_%d;',k)); %
    Hoop_Strain_EV(k,:,:) = Avg_Hoop_Strain;
 
%     eval(sprintf('load C:\\School_Work\\ESM_4714\\Long_Strain_Mode_%d;',k)); %
%     Long_Strain_EV(k,:,:) = Avg_Long_Strain;
 
end
 
clear Avg_Hoop_Strain
%clear Avg_Long_Strain
 
toc
 
%Node_Number = Hoop_Strain_EV(1,:,1);
Hoop_Strain_EV(:,:,1) = [];
 
HEV = Hoop_Strain_EV;
clear Hoop_Strain_EV;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Now we consider just the corner nodes 
%
% Midside nodes are elminiated
 
tic
 
load C:\\School_work\\ESM_4714\\All_Nodes_Elems_Position_in_Elems
 
NEP = All_Node_Elem_Pos;
clear All_Node_Elem_Pos
load_NEP = toc
 
counter = 0 
 
for n = 1:length(NEP)
    
    if (round(NEP(n,1,3)./2) - NEP(n,1,3)./2) > 0
       counter = counter+1;
       
        A(counter) = NEP(n,1,1);
        EV(:,counter) = HEV(:,n);        
                
    end
end
 
% node_Strain(1,:) = A;
% node_Strain(2:36,:) = EV;
 
EV_dim = size(EV);
 
% Clear extraneous variables
clear HEV
clear NEP
 
% set up freq vector in Hz, and in Rad.
 
    w = (2*pi).*[0.0029797 0.01862478 0.8229244 2.36246  13.50734 ... 
         14.45958  748.2332   748.285   1266.524 1266.526 ...
         1561.694  1759.436   1761.027  1761.04  1812.844 ...
         1813.000  1864.747   1865.151  2082.15  2293.602 ...
         2293.603  2454.926   2455.22   2554.576 2554.584 ...
         2575.635  2575.688   2673.731  2674.098 2916.426 ...
         2983.464  2984.044   3069.694  3069.71  3263.058];
 
freq = 0:0.5:2000;
w_samp = freq.*(2.*pi);
 
%@ w1 = 112 Hz spectral line 56
%@ w2 = 196 Hz spectral line 98
%@ w3 = 270 Hz spectral line 135
 
% Damping ratios calculated from real and imaginary components of FRF
 
d = 0.005*ones(1,length(w));
 
d( 7) = 0.00328;
d( 8) = 0.00328;
 
d( 9) = 0.003255;
d(10) = 0.003255;
 
d(13) = 0.00237;
d(14) = 0.00237;
 
d(17) = 0.00308;
d(18) = 0.00308;
 
 
%  Pack all the the variables into a temp directory
% cwd = pwd;
% cd C:\Temp;
% pack
% cd(cwd)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
no_modes = 18
 
tic
for j = 1:no_modes
    for k = 1:no_modes
        for u=1:1:length(w_samp)
            
%               Numerator(j,k,u) = (2.*(d(j)).*w(j)^2).*( transpose(EV(k,1:50))*EV(j,1:50));
%               Denominator(j,k,u) = ( (-w_samp(u).^2 + (w(j)).^2) + i.*(2.*(d(j)).*w(j)).*w_samp(u));
% 
%               rec_frf(j,k,u) = Numerator(j,k,u) / Denominator(j,k,u);
            rec_frf(:,:) = ((2.*(d(j)).*w(j)^2).*( EV(k,:).'*EV(j,:) ) / ( (-w_samp(u).^2 + (w(j)).^2) + i.*(2.*(d(j)).*w(j)).*w_samp(u)));
            
%            frf_3497(j,:,u) = rec_frf(852,:);
             frf_3497_temp(1,:,u) = rec_frf(852,:);
            
            if round (u./100) - u./100 == 0
                u
            end
            
        end
    end
    frf_3497 = squeeze(frf_3497_temp)
    clear rec_frf
                    
    j
    eval(sprintf('save C:\\School_work\\ESM_4714\\frf_3497_%d frf_3497',j));
    
end
toc 
 
%save C:\\School_work\\ESM_4714\\frf_3497_all frf_3497
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% These are the nodes that we are going to be looking at . . . 
 
c = 3;
 
e = 1;
f = 1;
g = 1;
h = 1;
 
m = 1;
n = 700;
o = 900;
p = 1634;
 
% Here are the plots
 
figure(1)
set(gcf, 'Position', [50 50 900 900])
subplot(2,1,1)
plot(freq(1:c), 20*log10(abs(squeeze(frf_3497(e,m,1:c)))),'-r')
hold on
plot(freq(1:c), 20*log10(abs(squeeze(frf_3497(f,m,1:c)))),'-g')
plot(freq(1:c), 20*log10(abs(squeeze(frf_3497(g,m,1:c)))),'-b')
plot(freq(1:c), 20*log10(abs(squeeze(frf_3497(h,m,1:c)))),'-k')
 
plot(freq(1:c), 20*log10(abs(squeeze(frf_3497(e,m,1:c) + frf_3497(f,m,1:c) + frf_3497(g,m,1:c) + frf_3497(h,m,1:c)))),'-om')
 
ylabel('Mag, dB')
title('1st column')
 
subplot(2,1,2)
plot(freq(1:c), (180/pi)*angle(squeeze(frf_3497(e,m,1:c))),'-r')
hold on
plot(freq(1:c), (180/pi)*angle(squeeze(frf_3497(f,m,1:c))),'-g')
plot(freq(1:c), (180/pi)*angle(squeeze(frf_3497(g,m,1:c))),'-b')
plot(freq(1:c), (180/pi)*angle(squeeze(frf_3497(h,m,1:c))),'-k')
 
plot(freq(1:c), (180/pi)*angle(squeeze(frf_3497(e,m,1:c) + frf_3497(f,m,1:c) + frf_3497(g,m,1:c) + frf_3497(h,m,1:c))),'-om')
ylabel('Phase, deg')
xlabel('Frequency, Hz')
legend('1,1','1,2','1,3','1,4','Mag 1')
 
 
 
% IFFT section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
total_resp(1,:) = squeeze(frf_3497(e,m,1:c) + frf_3497(f,m,1:c) + frf_3497(g,m,1:c) + frf_3497(h,m,1:c)).';
total_resp(2,:) = squeeze(frf_3497(e,n,1:c) + frf_3497(f,n,1:c) + frf_3497(g,n,1:c) + frf_3497(h,n,1:c)).';
total_resp(3,:) = squeeze(frf_3497(e,o,1:c) + frf_3497(f,o,1:c) + frf_3497(g,o,1:c) + frf_3497(h,o,1:c)).';
total_resp(4,:) = squeeze(frf_3497(e,p,1:c) + frf_3497(f,p,1:c) + frf_3497(g,p,1:c) + frf_3497(h,p,1:c)).';
 
% total_resp(5,:) = squeeze(frf_3497(1,1,1:c) + frf_3497(1,2,1:c) + frf_3497(1,3,1:c) + frf_3497(1,4,1:c)).';
 
freq = freq(1:c).';
 
figure(3)
set(gcf, 'Position', [50 50 900 900])
suptitle('5 node Fixed-Free Beam Frequency Response')
subplot(2,2,1)
plot(f(1:c), 20*log10(abs(total_resp(1,:))),'-om','Markersize',3)
ylabel('Mag of disp, dB')
title('node 2')
 
subplot(2,2,2)
plot(f(1:c), 20*log10(abs(total_resp(2,:))),'-om','Markersize',3)
ylabel('Mag of disp, dB')
title('node 3')
 
subplot(2,2,3)
plot(f(1:c), 20*log10(abs(total_resp(3,:))),'-om','Markersize',3)
ylabel('Mag of disp, dB')
xlabel('Freq, Hz')
title('node 4')
 
subplot(2,2,4)
plot(f(1:c), 20*log10(abs(total_resp(4,:))),'-om','Markersize',3)
ylabel('Mag of disp, dB')
xlabel('Freq, Hz')
title('node 5')
 
% figure(100)
% plot(f(1:c), 20*log10(abs(total_resp(5,:))),'-om','Markersize',3)
% ylabel('Mag of disp, dB')
% xlabel('Freq, Hz')
% title('IGNORE')
% 
 
dims = size(total_resp);
response = zeros(dims(1), 2*dims(2));
 
[response(1,:), time] = IFFT_time_resp(total_resp(1,1:c),freq(1:c));
[response(2,:), time] = IFFT_time_resp(total_resp(2,1:c),freq(1:c));
[response(3,:), time] = IFFT_time_resp(total_resp(3,1:c),freq(1:c));
[response(4,:), time] = IFFT_time_resp(total_resp(4,1:c),freq(1:c));
% [response(5,:), time] = IFFT_time_resp(total_resp(5,1:c),f(1:c));
 
 
figure(4)
set(gcf, 'Position', [50 50 900 900])
suptitle('5 node Fixed-Free Beam Time Response')
subplot(2,2,1)
plot(time, response(1,:),'-om','Markersize',3)
ylabel('disp, m')
title('node 2')
 
subplot(2,2,2)
plot(time, response(2,:),'-om','Markersize',3)
ylabel('disp, m')
title('node 3')
 
subplot(2,2,3)
plot(time, response(3,:),'-om','Markersize',3)
ylabel('disp, m')
xlabel('time, s')
title('node 4')
 
subplot(2,2,4)
plot(time, response(4,:),'-om','Markersize',3)
ylabel('disp, m')
xlabel('time, s')
title('node 5')
 
 
% figure(101)
% plot(time, response(5,:),'-om','Markersize',3)
% ylabel('disp, m')
% xlabel('time, s')
% title('IGNORE')