clc
clear all
close all

constants
weight=0;
j=0;
nmode=4;
%%
[phi E A]=findphi(weight,nmode);

for k=1:-0.005:0.0001
    k=k.^2;
    
    j=j+1;
    [B C E A]=(findBC(k,weight,nmode));
    
    F=A+1i*rho*c/2/k.*B+rho*(c/2/k)^2.*C;
    [G,~]=eig(E^-1*F);
    L=eig(E^-1*F);
    w1(:,j)=1./(sqrt(real(L)));
    gi(:,j)=imag(L)./real(L);
    v(:,j)=real(w1(:,j)*c/(2*k));
    
end
figure(1);
plot(v(1,:),gi(1,:),'*',v(2,:),gi(2,:),'*',v(3,:),gi(3,:),'*',v(4,:),gi(4,:),'*')
xlim([0 50])
ylim([-10 3])
grid on
legend ('1^s^t mode','2^n^d mode','3^r^d mode','4^t^h mode')
title('V g plot')
xlabel ('V [m/s]')
ylabel('g')

figure(2);
plot(v(1,:),w1(1,:)/2/pi,'*',v(2,:),w1(2,:)/2/pi,'*',v(3,:),w1(3,:)/2/pi,'*',v(4,:),w1(4,:)/2/pi,'*');
xlim([0 50])
grid on
legend ('1^s^t mode','2^n^d mode','3^r^d mode','4^t^h mode')
title('V \omega plot')
xlabel ('V [m/s]')
ylabel('\omega [Hz]')
