disp(' ')
disp(' dp_impedance.m   ver 1.1   October 10, 2006');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ')
%
clear f;
clear Z;
clear ff;
clear ZZ;
%
disp(' Enter structure type ');
disp(' 1=Beam longitudinal ')
disp(' 2=Rectangular Beam bending ')
disp(' 3=Square Tube Beam bending ')
disp(' 4=Plate edge ')
type=input(' ');
%
disp(' Enter rho (lbm/in^3)');
rho=input(' ');
rho=rho/386;  % change to lbf sec^2/in^4 
%
disp(' Enter CL (in/sec)');
CL=input(' ');
%
if(type==1 | type==2 )
   disp(' Enter cross-section area (in^2)');
   S=input(' ');
end
%
if(type==1)
    Z=rho*S*CL;
end
if(type==2)
    disp(' Enter thickness (in)');
    h=input(' ');
    for(i=10:10000)
        f(i)=i;
        Z(i)=0.67*rho*S*sqrt(CL*h*f(i))*sqrt(2);
    end
end
if(type==3)
    disp(' Enter width (in) ');
    w=input(' '); 
    disp(' Enter web thickness (in)');
    h=input(' ');
    S=4*w*h;
    E=(CL^2*rho);
    I=(1/12)*(w^4-(w-h)^4);
    for(i=10:10000)
        f(i)=i;
        omega=2*pi*i;
        Z(i)=sqrt(omega/2)*((rho*S)^(0.75))*((E*I)^(0.25));
    end
end
if(type==4)
    disp(' Enter thickness (in)');
    h=input(' ');   
    Z=CL*rho*h^2;
end
%
if(type==1 | type==4)
    disp(' ');
    out1=sprintf(' Z=%8.4g  lbf sec/in',Z);  
    disp(out1);
end
%
k=1;
mmm=size(f);
mn=max(mmm);
if(type==2 | type==3 ) 
    Q2=dp_plot(1,f,Z,'Frequency (Hz)','Impedance (lbf sec/in)','Driving Point Impedance  Beam in Bending Excited at Edge');    
    for(i=1:mn)
        if(Z(i)>0)
            ff(k)=f(i);
            ZZ(k)=1/Z(i);
            k=k+1;
        end
    end
end