%
disp(' ')
disp(' SEA_cylinder.m  ver 1.6  September 5, 2008')
disp(' By Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
%
%    m is the axial mode number
%    n is the circumferential mode number
%
clear f_bending;
clear k_bending;
%
clear nt;
clear ff;
clear fn;
clear ffb;
clear ma;
clear md;
clear mm;
clear nn;
clear psd;
clear a;
clear S;
clear SPL;
clear rho;
clear E;
clear K;
clear D;
clear h;
clear rem;
clear kcomp;
clear summary;
clear mtype;
%
nmax=170;
mmax=170;
%
disp(' ');
disp(' Select Type:  1=solid  2=honeycomb' );
itype=input(' ');
%
disp(' ');
if(itype==1)
    disp(' Enter material ' );
else
    disp(' Enter face sheet material ' );
end
%
disp('  1=aluminum  2=epoxy/graphite composite  3=other' );
%
imat=input(' ');
%
if(imat==1)
    E=9.9e+06;      % elastic modulus (psi)
    rho=0.0002539;  % mass/volume (lbf sec^2/in^4)
end
if(imat==2)
    E=10.0e+06;      % elastic modulus (psi)
    rho=0.00015024;  % mass/volume (lbf sec^2/in^4)       
end
if(imat > 2)
    disp(' ');
    disp(' Enter Elastic Modulus (lbf/in^2)' );
    E=input(' ');
    disp(' ');
    disp(' Enter mass/volume (lbm/in^3)' );
    rho=input(' ');
    rho=rho/386;
end
v=0.33;         % Poisson ratio
%
% damp=0.05;
disp(' ');
disp(' Enter damping ratio');
damp=input(' ');
loss=2*damp;
%
bc=0;
%
while(1)
    disp(' ');
    disp(' Selected Boundary Conditions ')
    disp(' 1=free-free ')   
    disp(' 2=fixed-fixed ')
    disp(' 3=fixed-free')
    disp(' 4=simply supported-simply supported ')
    bc=input(' ');
    if(bc==1 | bc==2 | bc==3 | bc==4)
        break;
    end
end
%
disp(' ');
disp(' Enter diameter (in)' );
diameter=input(' ');
a=diameter/2;
%
disp(' ');
disp(' Enter length (in)' );
L=input(' ');
%
disp(' ');
if(itype==1)
    disp(' Enter thickness (in)' );
    h=input(' ');
else
    disp(' Enter upper face sheet thickness (in)' );
    t1=input(' '); 
    disp(' Enter lower face sheet thickness (in)' );
    t2=input(' ');
    disp(' Enter the core thickness (in)');
    core=input(' ');
end
%
disp(' ');
disp(' Apply mass law at high frequencies?  1=yes 2=no ')
imlaw=input(' ');
%
length=L;
radius=a;
%
S=(pi*diameter*L);
%
if(itype==1)
    K=E*h/(1-v^2);
    D=E*h^3/(12*(1-v^2));
else
    K=E*(t1+t2)/(1-v^2); 
    h=core+(t1+t2)/2;
    D=E/(1-v^2);
    D=D*t1*t2*h^2/(t1+t2);
    h=(t1+t2);
end
%
DoverK=D/K;
%
out1=sprintf(' K=%8.4g  D=%8.4g  DoverK=%8.4g \n',K,D,DoverK);
disp(out1);
%
air_rho=1.2103e-007;   % lbf sec^2/in^4  
air_rho=2.8267e-008;   % at altitude
air_c=13500;           % in/sec
ci_air=air_rho*air_c;  % lbf sec/in^3
%
        CL=sqrt(E/rho);
    ring_f=CL/(pi*diameter);
clear arg;
if(itype==1)
    arg=h^3*rho/D;
else
    arg=h^3*rho/D;
end
% critical_f=(air_c^2/(2*pi*h))*sqrt(12*rho*(1-v^2)/E);
 critical_f=(air_c^2/(2*pi*h))*sqrt(arg);
%
%
mass=rho*S*h;
%
ma=zeros([3 3]);
ma(1,1)=1;
ma(2,2)=1;
ma(3,3)=1; 
%
ma=ma*(rho*h);
%
ma
out1=sprintf('\n mass=%8.4g  \n',mass);
disp(out1);
%
highest_AS=0;
%
kv=1;
%
progressbar % Create figure and set starting time 
disp(' ');
disp(' Calculating Modal Frequencies... ');
mstart=1;
if(bc==1)
    mstart=0;
end
for(m=mstart:mmax) % axial
%
    progressbar(m/mmax) % Update figure      
    if(bc==1 | bc==2)
        if(m==1)
            xo=[4.7 4.8];
            qcc=inline('cos(x).*cosh(x)-1','x');
            options=optimset('display','off');
            kx=fzero(qcc,xo,options);
        end
        if(m>=2)
            kx=(2*m+1)*pi/2;
        end
    end
    if(bc==3)
        if(m==1)
            xo=[1.8 1.9];
            qcc=inline('cos(x).*cosh(x)+1','x');
            options=optimset('display','off');
            kx=fzero(qcc,xo,options);
        end
        if(m>=2)
            kx=(2*m-1)*pi/2;
        end
    end
    if(bc==4)
        kx=m*pi;
    end
    if(m==0)
        kx=0.;
    end
    kx=kx/L;
    kx2=kx^2;
%    
    for(n=0:nmax) % circumferential
%
        clear ka;
%
        kt=n/a;
        kt2=kt^2;
        khat=sqrt(kx2+kt2);
%
        ka(1,1)=kx2+((1-v)/2)*kt2;
        ka(1,2)=(kt*kx*(1+v))/2; 
        ka(1,3)=(v/a)*kx; 
        ka(2,1)=ka(1,2);
        ka(2,2)=kt2+((1-v)/2)*kx2; 
        ka(2,3)=kt/a;
        ka(3,1)=ka(1,3);
        ka(3,2)=ka(2,3); 
        ka(3,3)=((DoverK)*(khat^4))+(1/a^2);    
        ka=K*ka;  
%
%  Optional for mode shapes
%
%        ka(1,2)=ka(1,2)*i;
%        ka(1,3)=ka(1,3)*i;
%        ka(2,1)=-ka(2,1)*i;
%        ka(3,1)=-ka(3,1)*i;
%
        [ModeShapes,som2]=eig(ka,ma);
        fn(1)=sqrt(som2(1,1));
        fn(2)=sqrt(som2(2,2));
        fn(3)=sqrt(som2(3,3));   
        fn=fn/(2.*pi);
        ff(kv)=fn(1);
        mm(kv)=m;
        nn(kv)=n;
%
        kn=n*pi/(pi*diameter);
        km=m*pi/length;
%
        if(m==0 & n==0)
            kn=1.0e-10;
            km=1.0e-10;
        end
%
        k=2*pi*fn(1)/air_c;
        if(k^2>=kn^2+km^2)
            nt(kv)=1;
            kratio=(km^2+kn^2)/k^2;
            rem(kv)=1/sqrt(1-kratio);
            
        else
            nt(kv)=0;
            rem(kv)=0;
        end
        if(m==0 & n==0)
            ff(kv)= ring_f;
            rem(kv)=1;
        end
%
        kv=kv+1;
%
%        out1=sprintf(' n=%d m=%d  fn1=%g  fn2=%g  fn3=%g',n,m,fn(1),fn(2),fn(3));
%        disp(out1)
    end
end
%
disp(' ');
disp(' Sorting... ');
clear fnm;
clear ntb;
clear nnb;
clear mmb;
clear remb;
clear re;
%
fnm=[ff',nn',mm',nt',rem'];
fnm=sortrows(fnm,1);
%
ffb=fnm(:,1);
nnb=fnm(:,2);
mmb=fnm(:,3);
ntb=fnm(:,4);
remb=fnm(:,5);
%
reff=zeros([1 kv]);
%
disp(' ')
disp('  fn(Hz)   n    m    kn    km   k   re  type  cb');
%
mmm=kv-1;
%
j=1;
%
progressbar; 
for(i=1:mmm)
    progressbar(i/mmm);
    if(ntb(i)==1)
       mtype='AF';     
    else
       mtype='AS';
       if(ffb(i)>highest_AS & (mmb(i)+nnb(i))~=0)
           highest_AS=ffb(i);
       end
    end
    kn=nnb(i)*pi/(pi*diameter);
    km=mmb(i)*pi/length;
    k=2*pi*ffb(i)/air_c;
    if(k>sqrt(kn^2+km^2))
        reff(i)=sqrt(1-(kn^2+km^2)/k^2);
    else
        reff(i)=0.01;
    end
%
    kcomp(i)=sqrt(kn^2+km^2);
%
    cw=0.;
    if(kcomp(i)>1.0e-10)
        cw=2*pi*ffb(i)/kcomp(i);
    end
    if(i<250)
        out1=sprintf('%d  %8.2f  %d  %d  %6.2g  %6.2g %6.2g %6.2g %s %10.5g',i,ffb(i),nnb(i),mmb(i),kn,km,k,reff(i),mtype,cw);
        disp(out1)
    end
%    
    if(i<=1000)
        summary(i,1)=i;
        summary(i,2)=ffb(i);
        summary(i,3)=nnb(i); 
        summary(i,4)=mmb(i);
        summary(i,5)=kn; 
        summary(i,6)=km; 
        summary(i,7)=kcomp(i);
        summary(i,8)=k;  
        summary(i,9)=reff(i);
        summary(i,10)=cw;
    end
%
    if(nnb(i)==0)
       f_bending(j)=ffb(i);
       k_bending(j)=kcomp(i);
       j=j+1;
    end
%
    beta=sqrt(h^2/(12*a^2));
    ksab=kn*a*sqrt(beta);
    kzab=km*a*sqrt(beta);
%    
    yup=0;
    ylow=0;
    if(ksab<=0.5)
         yup=0.5+sqrt(0.5^2-ksab^2);
         ylow=0.5-sqrt(0.5^2-ksab^2);
    end
%  
    if((kzab< ylow | kzab> yup) | ksab>0.5)
%        f_bending(j)=ffb(i);
%        k_bending(j)=kcomp(i);
%        j=j+1; 
    end
%
end
%
disp(' ')
disp(' Candidate Critical Frequencies ')
progressbar;
for(i=1:mmm)
    progressbar(i/mmm);
    if(kcomp(i)>1.0e-10)
        cw=2*pi*ffb(i)/kcomp(i);   
        kn=nnb(i)*pi/(pi*diameter);
        km=mmb(i)*pi/length;   
        k=2*pi*ffb(i)/air_c;   
        if( abs(k-kcomp(i))<0.05*k)
            out1=sprintf('%d  %8.2f  %d  %d  %6.2g  %6.2g %6.2g %6.2g %s %10.5g',i,ffb(i),nnb(i),mmb(i),kn,km,k,reff(i),mtype,cw);
            disp(out1)  
        end
    end
end
%
out1=sprintf('\n ring_f=%10.5g Hz   flat plate critical_g=%10.5g Hz ',ring_f,critical_f);
disp(out1)
%
out1=sprintf('\n Highest AS frequency = %10.5g Hz \n',highest_AS);
disp(out1);  
%
clear fc;
clear nn;
%
ioct=2;
if(ioct==2)
	fc(1)=2.5;
	fc(2)=3.15;
	fc(3)=4.;;
	fc(4)=5.;
	fc(5)=6.3;
	fc(6)=8.;
	fc(7)=10.;
	fc(8)=12.5;
	fc(9)=16.;
	fc(10)=20.;
	fc(11)=25.;
	fc(12)=31.5;
	fc(13)=40.;
	fc(14)=50.;
	fc(15)=63.;
	fc(16)=80.;
	fc(17)=100.;
	fc(18)=125.;
	fc(19)=160.;
	fc(20)=200.;
	fc(21)=250.;
	fc(22)=315.;
	fc(23)=400.;
	fc(24)=500.;
	fc(25)=630.;
	fc(26)=800.;
	fc(27)=1000.;
	fc(28)=1250.;
	fc(29)=1600.;
	fc(30)=2000.;
	fc(31)=2500.;
	fc(32)=3150.;
	fc(33)=4000.;
	fc(34)=5000.;
	fc(35)=6300.;
	fc(36)=8000.;
	fc(37)=10000.;
	fc(38)=12500.;
	fc(39)=16000.;
	fc(40)=20000.;
    imax=40;
end
%
nun=zeros([1 imax]);
rem_ave=zeros([1 imax]);
rem_sum=zeros([1 imax]);
nrband=zeros([1 imax]);
nslow=zeros([1 imax]);
nfast=zeros([1 imax]);
%
ocd=2^(1/6);
for(j=1:imax)
    fl=fc(j)/ocd;
    fu=fc(j)*ocd;
    for(i=1:(kv-1))
%
        if(ffb(i)>=fl & ffb(i) < fu )   
           if(ntb(i)==1)
               nun(j)=nun(j)+1;
               nfast(j)=nfast(j)+1;
           else
               nun(j)=nun(j)+1;
               nslow(j)=nslow(j)+1;   
           end
           rem_sum(j)=rem_sum(j)+reff(i);
        end
        if(ffb(i)>fu)
            break;
        end
    end
    rem_ave(j)=0;
    if(nun(j)>0)
       rem_ave(j)=rem_sum(j)/nun(j);
    end
    nun(j)=nun(j)/(fu-fl);
end  
%
md=[fc' nun'];
md=sortrows(md,1);
%
Q2=SEA_plot(1,md(:,1),md(:,2),'Frequency (Hz)','Modes/Hz','Cylinder Modal Density',1);
%
disp(' ');
disp(' Write modal density to ASCII file? ')
pchoice=input(' 1=yes 2=no ');
if(pchoice == 1) 
%
      clear writedata;
      [writefname, writepname] = uiputfile('*','Save modal density data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [ md(:,1), md(:,2) ];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %g  %g \n',writedata');
	  fclose(fid);
end  
%
disp(' ');
disp('  Power Transmissibility Analysis ');
%
%   Power Transmissibility Analysis
%
clear trans;
clear ptrans;
clear pfc;
clear Rratio;
clear re;
%
k=1;
for(i=1:imax)
%    
    omega=2*pi*fc(i);
    if(nun(i)>0)
%        Z=4*mass/nun(i);
%        Z=2*mass*sqrt(omega*loss/nun(i))     
%        Z=2.3*CL*rho*h^2;
        Z=rho*h*omega*loss;
    else
        Z=1.0e+80;
    end
%
%    v0=fc(i)/ring_f;
%    radiation_efficiency=0.346*sqrt(1-v^2)*v0*ring_f/critical_f;
%    
    radiation_efficiency=rem_ave(i);
%    
%    if(radiation_efficiency>1)
%        radiation_efficiency=1;
%    end
%
    Rmech=Z;
    Rrad=ci_air*radiation_efficiency;
%
    RRR=Rrad/(Rrad+Rmech);    
%  
%    out1=sprintf(' fc=%8.3f  RRR=%8.3g  Rrad=%8.4g  Rmech=%8.4g  Reff=%8.4g\n',fc(i),RRR,Rrad,Rmech,radiation_efficiency);    
%    disp(out1)    
%
    alpha=pi;
%
    trans(i)=alpha*nun(i)*air_c^2/(ci_air*mass);
    trans(i)=trans(i)*RRR;
%
    if(imlaw==1)
%
%   Mass Law + 4.5 dB
%
        ML=(S/mass);   % Mass law
        ML=ML^2;
        ML=ML*2.82;
        if(fc(i)>1.2*critical_f | fc(i)>2.*ring_f)
            trans(i)=ML;
        end
 %   
    end
%    
    trans(i)=trans(i)/386^2;         
%
%   Smooth
%
    if(i>=3 & (trans(i-1)-trans(i-2))>0 )
        if(trans(i)<trans(i-1)*0.92 )
            trans(i)=trans(i-1)*0.92; 
        end
    end
    if(i>1 & trans(i)<trans(i-1)*0.7 ) 
        trans(i)=trans(i-1)*0.7;
    end  
    if(trans(i)>1e-08 & fc(i)>10)
        pfc(k)=fc(i);
        Rratio(k)=RRR;
        ptrans(k)=trans(i);
        re(k)=radiation_efficiency; 
        k=k+1;
    end
%
%   out1=sprintf(' fc=%8.3f  RRR=%g   \n',fc(i),RRR);    
%   disp(out1)
%
end
%
Q2=SEA_plot(2,pfc,re,'Frequency (Hz)',' ','Radiation Efficiency',1);
Q2=SEA_plot(3,pfc,Rratio,'Frequency (Hz)',' ','Impedance Ratio',1);
Q2=SEA_plot(4,pfc,ptrans,'Frequency (Hz)','G^2 / psi^2','Power Transmissibility',1);
Q2=SEA_plot(5,kcomp,ffb,'Wavenumer (rad/in)','Frequency (Hz)','Frequency vs. Wavenumber',3);
%
disp(' ');
disp(' Write Transmissibility (G^2/psi^2) to ASCII file? ')
pchoice=input(' 1=yes 2=no ');
if(pchoice == 1) 
%
      clear writedata;
      [writefname, writepname] = uiputfile('*','Save Transmissibility data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [ pfc',  ptrans' ];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %g  %g \n',writedata');
	  fclose(fid);
end 
%
%
disp(' ');
disp(' Write Frequency vs. Wavenumber to ASCII file? ')
pchoice=input(' 1=yes 2=no ');
if(pchoice == 1) 
%
      clear writedata;
      [writefname, writepname] = uiputfile('*','Frequency vs. Wavenumber');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [ kcomp',  ffb ];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %g  %g \n',writedata');
	  fclose(fid);
end 
%
%
disp(' ');
disp(' Write Full Text Ouput to ASCII file? ')
pchoice=input(' 1=yes 2=no ');
%
ss=size(summary);
if(ss(2)>ss(1))
    summary=summary';
end
%
%
if(pchoice == 1) 
%
      clear writedata;
      [writefname, writepname] = uiputfile('*','Full Text Output');
	  writepfname = fullfile(writepname, writefname);
%	  writedata = [ summary,  mtype ];
	  fid = fopen(writepfname,'w');
%
      ss=size(summary);
      for(i=1:ss(1)) 
%
        if(summary(i,7)>=summary(i,8))
            mtype='AS';     
        else
            mtype='AF';
        end
%
        fprintf(fid,' %d \t %8.4g \t %d \t %d \t %8.4g \t %8.4g \t %8.4g \t %8.4g \t %8.4g \t %10.5g \t %s\n',summary(i,1),summary(i,2),summary(i,3),summary(i,4),summary(i,5),summary(i,6),summary(i,7),summary(i,8),summary(i,9),summary(i,10),mtype );
%        
      end
%     
	  fclose(fid);
end 
%
disp(' ');
disp(' Perform force response analysis?  1=yes  2=no ');
nop=input(' ')
if(nop==1)
disp('  Forced Response Analysis ');
%
%    Forced Response Analysis
%
disp(' ');
disp(' Enter one-third octave SPL(dB) input file.')
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
file_choice = input('');
%
clear SPL;
%
if(file_choice==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
%      
    fid = fopen(filename,'r');
    SPL = fscanf(fid,'%g %g',[2 inf]);
    SPL=SPL';
else
    SPL = input(' Enter the matrix name:  ');
end
%
n=size(SPL);
%
reference=1;
ms=0;
%
for(i=1:n(1))
    pressure_ms=reference*(10.^(SPL(i,2)/10.) );
    ms=ms+pressure_ms;
end    
%
dbrms=10*log10(ms/reference)
%
out1=sprintf('SOUND PRESSURE LEVEL   OASPL = %8.3g dB  (Ref 20 micro Pa)',dbrms);
%
Q2=SEA_plot(6,SPL(:,1),SPL(:,2),'Center Frequency (Hz)','SPL (dB)',out1,2);
%
clear pressure_rms;
clear pressure_psd;
clear psd;
clear fff;
clear NL;
%
reference=2.9e-09;
%
NL=size(SPL(:,1));
%
nnmm=NL(1);
%
for(i=1:nnmm)
     pressure_rms(i)=reference*(10.^(SPL(i,2)/20.) );
     ffc=SPL(i,1);
     ffu=ffc*2^(1/6);
     ffl=ffc/2^(1/6);
     ddf(i)=ffu-ffl;
     pressure_psd(i)=(pressure_rms(i))^2/ddf(i);
end    
%
clear psd;
clear fff;
clear ipmax;
k=1;
ipmax=size(ptrans);
for(i=1:ipmax(2))
    for(j=1:NL-1)
        if(pfc(i)>=SPL(j,1) & pfc(i)<SPL((j+1),1) )
%          
            psd(k)=ptrans(i)*pressure_psd(j);          
            fff(k)=pfc(i);
            k=k+1;
        end
    end
end
%
Q2=SEA_plot(7,fff,psd,'Frequency (Hz)','Accel (G^2/Hz)','Power Spectral Density',1);
%
clear fq;
clear aq;
k=1;
%
clear nnpp;
nnpp=size(psd);
if(nnpp<=1)
    disp(' nnpp error ');
end
for(i=1:nnpp(2))
    if(psd(i)>1.0e-08 & fff(i)<=2000)
        fq(k)=fff(i);
        aq(k)=psd(i);
        k=k+1;
    end
end
%
[s,grms] = calculate_PSD_slopes(fq,aq);
%
out5 = sprintf('\n Overall Acceleration = %10.3g GRMS up to 2000 Hz',grms);
disp(out5);
%
disp(' ');
disp(' Write response PSD to ASCII file? ')
pchoice=input(' 1=yes 2=no ');
if(pchoice == 1) 
%
      clear writedata;
      [writefname, writepname] = uiputfile('*','Save PSD data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [ fff',  psd' ];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %g  %g \n',writedata');
	  fclose(fid);
end  
end