%
disp(' ');
disp(' ss_panel_flutter.m  ver 1.0   August 24, 2010 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' Movchan (157) formula for a panel simply supported on all edges.');
disp(' ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=other ');
%
if(imat==1)
	E=1.0e+07;
    mu=0.3;
end
if(imat==2)
    E=3.0e+007;
    mu=0.3;
end
if(imat ~=1 && imat ~=2)
    disp(' ');   
    E = input(' Enter the elastic modulus (lbf/in^2) ');
    disp(' ');
    mu = input(' Enter the Poisson ratio ');           
 end
%
Lx=input(' Enter length (in) ');
Ly=input(' Enter width (in) ');
%
L_ratio=Lx^2/Ly^2;
%
h=input(' Enter thickness (in) ');
%
D=E*h^3/(12*(1-mu^2));
%
while(1)
%
    alt_ft=input(' Enter the altitude (feet) ');
    alt_m=alt_ft*0.0003048;
    if(alt_m<=20)
        break;
    else
        disp(' Maximum = 65000 ft ' );
    end
end    
%
alt(1)=	0;	md(1)=	1.226;	cc(1)=	340.2;
alt(2)=	1;	md(2)=	1.112;	cc(2)=	336.3;
alt(3)=	2;	md(3)=	1.007;	cc(3)=	332.4;
alt(4)=	3;	md(4)=	0.9096;	cc(4)=	328.5;
alt(5)=	4;	md(5)=	0.8195;	cc(5)=	324.5;
alt(6)=	5;	md(6)=	0.7365;	cc(6)=	320.4;
alt(7)=	6;	md(7)=	0.66;	cc(7)=	316.3;
alt(8)=	7;	md(8)=	0.5898;	cc(8)=	312.1;
alt(9)=	8;	md(9)=	0.5254;	cc(9)=	307.9;
alt(10)=	9;	md(10)=	0.4666;	cc(10)=	303.7;
alt(11)=	10;	md(11)=	0.4129;	cc(11)=	299.3;
alt(12)=	11;	md(12)=	0.3641;	cc(12)=	295;
alt(13)=	12;	md(13)=	0.3104;	cc(13)=	295;
alt(14)=	13;	md(14)=	0.2652;	cc(14)=	295;
alt(15)=	14;	md(15)=	0.2266;	cc(15)=	295;
alt(16)=	15;	md(16)=	0.1936;	cc(16)=	295;
alt(17)=	16;	md(17)=	0.1654;	cc(17)=	295;
alt(18)=	17;	md(18)=	0.1413;	cc(18)=	295;
alt(19)=	18;	md(19)=	0.1207;	cc(19)=	295;
alt(20)=	19;	md(20)=	0.1032;	cc(20)=	295;
alt(21)=	20;	md(21)=	0.0881;	cc(21)=	295;
%
for(i=1:20)
%    
    if(alt_m==alt(i))
        mdens=md(i);
        c=cc(i);
        break;
    end
%
    if(alt_m==alt(i+1))
        mdens=md(i+1);
        c=cc(i+1);
        break;
    end
%
    if(alt_m>alt(i) && alt_m<alt(i+1))
        length=alt(i+1)-alt(i);
            xx=alt_m-alt(i);
            c2=xx/length;
            c1=1-c2;
         mdens=c1*md(i)+c2*md(i+1);
             c=c1*cc(i)+c2*cc(i+1);        
    end
%
end
mdens=mdens*9.36e-008;
c=c*39.37;
rho_c2=mdens*c^2;
%
out1=sprintf('\n  Plate Stiffness Factor = %8.4g lbf in ',D);
out2=sprintf('\n  Speed of sound = %8.4g in/sec  ',c);
  out3=sprintf('                 = %8.4g ft/sec  ',c/12);
out4=sprintf('\n     Air density = %8.4g lbf sec^2/in^4 ',mdens);
  out5=sprintf('                 = %8.4g lbm/in^3      ',mdens*386.04);
out6=sprintf('\n          rho_c2 = %8.4g lbf/in^2 ',rho_c2);
%
disp(out1);
disp(out2);
disp(out3);
disp(out4);
disp(out5);
disp(out6);
%
num=D*8*pi^3;
den=rho_c2*(Lx^3)*3*sqrt(3);
term=(5+L_ratio)*sqrt(2+L_ratio);
Mcr  = (num/den)*term;
%
out1=sprintf('\n Mcr = %8.4g ',Mcr);
disp(out1);
%