%2 degree of freedom flutter solver
mass=0.0486;
polarI=4.05e-5;
wh=12.0
wh=wh*2.0*pi
wa=105.0
wa=wa*2.0*pi
Kh=wh*wh*mass
Ka=wa*wa*polarI
rho=1.225;
S=0.1*.45;
chord=0.1;
Cla=2.0*pi;
xac=0.25*chord;
xcg=-0.01*chord;
crossI=xcg*mass;
Vmax=input('Maximum Velocity = ');
%
%
V=0.0;
dv=Vmax/100.0;
for i=1:100,
 A=mass*polarI-crossI*crossI;
 B=0.5*rho*V*V*S*Cla*(crossI+xac*mass)-Ka*mass-Kh*polarI;
 C=Kh*(Ka-0.5*rho*V*V*S*Cla*xac);
%
%
 w1=(-B+sqrt(B*B-4.0*A*C))/(2.0*A);
 w2=(-B-sqrt(B*B-4.0*A*C))/(2.0*A);
 w11(i)=sqrt(w1)/2.0/pi;
 w22(i)=sqrt(w2)/2.0/pi;
 u(i)=V;
 V=V+dv;
end
w33=imag(w11);
w44=imag(w22);
plot(u,w11,'-w',u,w22,'-w');
xlabel('Velocity (m/s)');
ylabel('Frequency (rad/s)');
title('Flutter Prediction Plot');