% Exercise 6: Comparasion between Techniques of Fixed Points and 
%             Non Linear Optimizaton Techniques - Viscous Model 
clc
clear
%
% Global variables
global Plothandle m1 k1 c1 ksi1 m2 omegaot ome1 ome2 C11 omega iomegmax
%
%Primary system data
m1=0.191;
k1=17.4e4;
ksi1=0.01;
c1=ksi1*2*m1*sqrt(k1/m1);
mu=0.1;
m2=mu*m1;

% Band of frequency to optimal design (change to broadband or narrowband)
  ome1=100;
  ome2=2000;
%   ome1=930;
%   ome2=970;
%    ome1=1200;
%    ome2=1250;
C11=(ome1+ome2)/2;
%
omegaot=ome1:0.5:ome2;
%
% Band of frequency to build the graph (view results)
omega2=2500;
omega1=2*pi*0;
%
omega=omega1:0.5:omega2;
%
% FRF of the primary system (receptance)
for iomeg=1:length(omega)
    beta(iomeg)=omega(iomeg)/sqrt(k1/m1);
    H1(iomeg)=1./(-m1*(omega(iomeg))^2+k1+2*m1*ksi1*sqrt(k1/m1)*omega(iomeg)*i);
end
%
subplot(2,1,1)
plot(beta,abs(H1),'k-')
%axis tight
axis([0.5 1.5 0 5e-5])
xlabel('omega/omegan')
ylabel('abs(H1)')
title('FPT and NLOT, with ksi1=0')
%
pause   % use enter to continue
%
hold on
%
Plothandle=plot(beta,abs(H1),'r-','EraseMode','xor')
%
% Non Linear Optimization Technics - NLOT
lam=[1000 100];                             % initial point of search
lambda=fminsearch('fobj_exercise6',lam);   % otpimization
%
% Technics of Fixed POints - FPT
% The optimal parameters using FPT
alphao=1/(1+mu);
k2=m2*(alphao*sqrt(k1/m1))^2;
Xi2o=1/alphao*sqrt(3*mu/(8*(1+mu)^3));
c2=Xi2o*2*m2*sqrt(k2/m2);
%
 for iomeg=1:length(omega)
     H11(iomeg)=(-omega(iomeg)^2*m2+k2+i*omega(iomeg)*c2)/((-omega(iomeg)^2*m2+k2)*(-omega(iomeg)^2*m1+k1)-omega(iomeg)^2*m2*k2+i*omega(iomeg)*c2*(-omega(iomeg)^2*m1+k1-omega(iomeg)^2*m2));
 end
%
hold on
plot(beta,abs(H11),'k-')
%
% % Considering ksi2=0
c2=0;

 for iomeg=1:length(omega)
     H11(iomeg)=(-omega(iomeg)^2*m2+k2+i*omega(iomeg)*c2)/((-omega(iomeg)^2*m2+k2)*(-omega(iomeg)^2*m1+k1)-omega(iomeg)^2*m2*k2+i*omega(iomeg)*c2*(-omega(iomeg)^2*m1+k1-omega(iomeg)^2*m2));
 end
 
plot(beta,abs(H11),'b-')


% Considering ksi2=infinit
c2=1000000000000;

 for iomeg=1:length(omega)
     H11(iomeg)=(-omega(iomeg)^2*m2+k2+i*omega(iomeg)*c2)/((-omega(iomeg)^2*m2+k2)*(-omega(iomeg)^2*m1+k1)-omega(iomeg)^2*m2*k2+i*omega(iomeg)*c2*(-omega(iomeg)^2*m1+k1-omega(iomeg)^2*m2));
 end
 
plot(beta,abs(H11),'g-')



legend('without absorber','NLOT','FPT')
%
%
% Primary system with damping ksi1=0.05
%
ksi1=0.1;
c1=ksi1*2*m1*sqrt(k1/m1);
%
% FRF of the primary systema with damping
for iomeg=1:length(omega)
    beta(iomeg)=omega(iomeg)/sqrt(k1/m1);
    H1(iomeg)=1./(-m1*(omega(iomeg))^2+k1+2*m1*ksi1*sqrt(k1/m1)*omega(iomeg)*i);   %Receptância
end
%
subplot(2,1,2)
plot(beta,abs(H1),'k-')
axis([0.5 1.5 0 5e-5])
xlabel('omega/omegan')
ylabel('abs(H1)')
title('FPT and NLOT with damped primary system')
%
pause   % click enter to continue
%
hold on
%
Plothandle=plot(beta,abs(H1),'r-','EraseMode','xor')
%
% NLOT
lam=[150 100];
lambda2=fminsearch('fobj2_exercise6',lam);
%
% FPT
% Optimal parameters obtained considering ksi1=0
alphao=1/(1+mu);
k2=m2*(alphao*sqrt(k1/m1))^2;
Xi2o=1/alphao*sqrt(3*mu/(8*(1+mu)^3));
c2=Xi2o*2*m2*sqrt(k2/m2);
%
 for iomeg=1:length(omega)
       H11(iomeg)=(-omega(iomeg)^2*m2+k2+i*omega(iomeg)*c2)/((-omega(iomeg)^2*m2+k2)*(-omega(iomeg)^2*m1+k1+i*omega(iomeg)*c1)- ...
                  omega(iomeg)^2*m2*k2+i*omega(iomeg)*c2*(-omega(iomeg)^2*m1+k1+i*omega(iomeg)*c1-omega(iomeg)^2*m2));
 end
%
plot(beta,abs(H11),'b-')
legend('Without absorber','NLOT','FPT')
