% descontinuidades na FRF do oscilador de Duffing
clc
clearvars
close all

% parâmetros do sistema não linear cúbico (W&T, pp. 543-545)
m=1;
k=1;
k3=1;
c=0.3; % 0.01, 0.03, 0.1, 0.3

% dados da excitação
X=100; % 0.1 0.3 1.0 3.0 10.0 30.0 100.0 (W&T, pp. 543-545)

% pontos de bifurcação
  % coeficientes do polinômio de grau 10 (discriminante)
ca=256/(729*k3^6);
c10=-64*c^2*m^4;
c8=-128*c^4*m^2+256*c^2*k*m^3;
c6=-64*c^6+256*c^4*k*m-384*c^2*k^2*m^2+48*k3*m^3*X^2;
c4=-128*c^4*k^2+256*c^2*k^3*m+432*c^2*k3*m*X^2-144*k*k3*m^2*X^2;
c2=-64*c^2*k^4-432*c^2*k*k3*X^2+144*k^2*k3*m*X^2;
c0=-48*k^3*k3*X^2-243*k3^2*X^4;
  % montagem do vetor de coeficientes
cp1=ca*[c10 0 c8 0 c6 0 c4 0 c2 0 c0];
  % determinação das raízes
rp1=roots(cp1);
  % seleção dos pontos de interesse
pb=[];
for ir=1:length(rp1)
  if (imag(rp1(ir)) == 0 && real(rp1(ir)) > 0)
    pb=[pb rp1(ir)];
  end
end
  % visualização de raízes selecionadas (frequências de salto)
pb=sort(pb)

% montagem de FRF
  % faixa de frequência (escolha em função dos parâmetros do sistema)
om=1:0.25:30;
  % inicialização de amplitudes de resposta e ...?
Y1=[]; om1=[]; Y2=[]; om2=[]; Y3=[]; om3=[];
for im=1:length(om)
  % coeficientes do polinômio de grau 6 em Y (W&T, p. 83, eq. 3.10)
  cc6=0.75^2*k3^2;
  cc4=-1.5*m*k3*om(im)^2+1.5*k*k3;
  cc2=m^2*om(im)^4-2*m*k*om(im)^2+k^2+c^2*om(im)^2;
  cc0=-X^2;
  % montagem do vetor de coeficientes em Y
  cp2=[cc6 0 cc4 0 cc2 00 cc0];
  % determinação das raízes
  rp2=roots(cp2);
  % seleção das amplitudes Y de interesse
  Yaux=[];
  for ir=1:length(rp2)
    if (imag(rp2(ir)) == 0 && real(rp2(ir)) > 0)
      Yaux=[Yaux rp2(ir)];
    end
  end
  Yaux=sort(Yaux);
  % montagem dos vetores de amplitude
  if length(Yaux) == 1
    Y1=[Y1 Yaux];
    om1=[om1 om(im)];
  else
    Y1=[Y1 Yaux(3)];
    om1=[om1 om(im)];
    Y2=[Y2 Yaux(2)];
    om2=[om2 om(im)];
    Y3=[Y3 Yaux(1)];
    om3=[om3 om(im)];
  end
end
  % cômputo dos valores da FRF (W&T, p. 84, eq. 3.13)
Lb1=1./(k+0.75*k3*Y1.^2-m*om1.^2+i*c*om1);
Lb2=1./(k+0.75*k3*Y2.^2-m*om2.^2+i*c*om2);
Lb3=1./(k+0.75*k3*Y3.^2-m*om3.^2+i*c*om3);

% visualização de FRF (magnitude e fase)
  % magnitude
figure(1)
plot(om1,abs(Lb1),'.g')
hold on
plot(om2,abs(Lb2),'.r')
plot(om3,abs(Lb3),'.b')
xlabel('frequência (rad/s)')
ylabel('magnitude (m/N)')
grid on
  % fase
figure(2)
plot(om1,atan2(imag(Lb1),real(Lb1)),'.g')
hold on
plot(om2,atan2(imag(Lb2),real(Lb2)),'.r')
plot(om3,atan2(imag(Lb3),real(Lb3)),'.b')
xlabel('frequência (rad/s)')
ylabel('fase (rad)')
grid on
  % Nyquist
figure(3)
plot(Lb1,'.g')
hold on
plot(Lb2,'.r')
plot(Lb3,'.b')
axis equal
xlabel('Re(\Lambda)')
ylabel('Im(\Lambda)')
grid on