disp(' ')
disp(' vrs_force.m   ver 1.2   December 20, 2013')
disp(' by Tom Irvine   Email: tom@vibrationdata.com')
disp(' ')
disp(' This program calculates the response of a single-degree-of-freedom')
disp(' system to a force power spectral density where the natural frequency ')
disp(' is an independent variable. ')
disp(' ')
disp(' The power spectral density may vary arbitrarily with frequency. ')
disp(' ')
%
clear THM;
clear F;
clear a;
clear fi;
clear fn;
clear omega;
clear omega2;
clear omegan;
clear omegan2;
clear ai;
clear length;
%
close all;
%
disp(' ');
disp(' Select units:   1=English  2=metric ');
iu=input(' ');
%
fig_num=1;
%
if(iu==1)
    disp(' The input file must have two columns: frequency(Hz) & force(lbf^2/Hz)')
else
    disp(' The input file must have two columns: frequency(Hz) & force(N^2/Hz)')    
end
%
[THM]=enter_PSD();
%
tpi=2.*pi;
%
F=THM(:,1);
a=THM(:,2);
%
nnn=length(F);
%
f1=F(1);
f2=F(nnn);
%
[m0,frms]=cm0(a,F,nnn);
%
figure(fig_num);
fig_num=fig_num+1;
plot(F,a);
grid on;
xlabel('Frequency (Hz)');
if(iu==1)
    ylabel('Force (lbf^2/Hz)');
else
    ylabel('Force (N^2/Hz)');   
end    
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
if(iu==1)
    out1=sprintf('Force PSD  Overall Level = %10.3f lbf RMS',frms);
    out5 = sprintf('\n Input PSD Overall Level = %10.3f lbf RMS',frms);
else
    out1=sprintf('Force PSD  Overall Level = %10.3f N RMS',frms);
    out5 = sprintf('\n Input PSD Overall Level = %10.3f N RMS',frms);    
end    
%
%
title(out1);
disp(out5);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
if(iu==1)
    mass=input('  Enter mass (lbm)  ');
else
    mass=input('  Enter mass (kg)  ');    
end
%
mass_orig=mass;
if(iu==1)
    mass=mass/386;
end    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
Q=input(' Enter the amplification factor (typically Q=10) ');
damp=1/(2*Q);
%
disp(' ');
disp('  Also calculate fatigue damage spectrum?  1=yes 2=no  ');
ird=input(' ');
disp(' ');
%
if(ird==1)
    bex=input(' Enter fatigue exponent  ');
    duration=input(' Enter duration (sec)  ');
    disp(' ');
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ocn=1./24.;
%
fi(1)=f1;
%
j=2;
while(1)
    fi(j)=fi(j-1)*2^(ocn);
    if(fi(j)>f2)
        fi(j)=f2;
        break;
    end
    j=j+1;
end
fn=fi;
%
n_ref=length(fi);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[ai]=interp_psd_oct(F,a,fi);
%
%
omega=zeros(n_ref,1);
omega2=zeros(n_ref,1);
omegan=zeros(n_ref,1);
omegan2=zeros(n_ref,1);
stiffness=zeros(n_ref,1);
%
for j=1:n_ref
%    
    omega(j)=tpi*fi(j);
    omega2(j)=omega(j)^2;
%
    omegan(j)=tpi*fn(j);
    omegan2(j)=omegan(j)^2;
%
    stiffness(j)=omegan2(j)*mass;
%
end
%
drms=zeros(n_ref,1);
vrms=zeros(n_ref,1);
grms=zeros(n_ref,1);
FT_rms=zeros(n_ref,1);
a_damage=zeros(n_ref,1);
ft_damage=zeros(n_ref,1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
progressbar;
for ijk=1:n_ref
    progressbar(ijk/n_ref);
%
    dpsd=zeros(n_ref,1);
    vpsd=zeros(n_ref,1);
    apsd=zeros(n_ref,1);
    FT_psd=zeros(n_ref,1);
%
%
%      Transmitted Force
%
    for j=1:n_ref
%
        rho=fi(j)/fn(ijk);
        num= 1.+(2*damp*rho)^2;
        den= (1-rho^2)^2 + (2*damp*rho)^2;
        t= num/den;
        FT_psd(j)=( t*ai(j) );
%
    end
    [FT_m0,FT_rms_ind]=cm0(FT_psd,fi,n_ref);
    FT_rms(ijk)=FT_rms_ind;
%    
%      Absolute Displacement
%
    for j=1:n_ref
%        
        num=( (omegan2(ijk)/stiffness(ijk))^2 );
        den=( (omegan2(ijk)-omega2(j))^2) + ( (2*damp*omega(j)*omegan(ijk))^2);
        t= num/den;
        dpsd(j)=( t*ai(j) );
%
    end
    [m0,rms]=cm0(dpsd,fi,n_ref);
    drms(ijk)=rms;
%
%      Velocity
%
    for j=1:n_ref
%     
        num=( (omega(j)*omegan2(ijk)/stiffness(ijk))^2 );
        den=( (omegan2(ijk)-omega2(j))^2) + ( (2*damp*omega(j)*omegan(ijk))^2);
        t= num/den;
        vpsd(j)=( t*ai(j) );
%
    end
    [m0,rms]=cm0(vpsd,fi,n_ref);
    vrms(ijk)=rms;
%
%      Acceleration
%
    for j=1:n_ref
%
        num=( (omega2(j)*omegan2(ijk)/stiffness(ijk))^2 );
        den=( (omegan2(ijk)-omega2(j))^2) + ( (2*damp*omega(j)*omegan(ijk))^2);
        t= num/den;
        apsd(j)=( t*ai(j) );
%
    end
    if(iu==1)
        apsd=apsd/(386^2);
    else
        apsd=apsd/(9.81^2);        
    end 
    [m0,rms]=cm0(apsd,fi,n_ref);
    grms(ijk)=rms;
%
%
    if(ird==1)
        [m1]=cm1(apsd,fi,n_ref);
        [m2]=cm2(apsd,fi,n_ref);
        [m4]=cm4(apsd,fi,n_ref);
%
        EP=sqrt(m4/m2);
%
        rms=grms(ijk);
        [aqd]=Dirlik_basic(duration,bex,m0,m1,m2,m4,EP,rms);
        a_damage(ijk)=aqd;
%
        [FT_m1]=cm1(FT_psd,fi,n_ref);
        [FT_m2]=cm2(FT_psd,fi,n_ref);
        [FT_m4]=cm4(FT_psd,fi,n_ref);
%
        FT_EP=sqrt(FT_m4/FT_m2);
%
        [fqd]=Dirlik_basic(duration,bex,FT_m0,FT_m1,FT_m2,FT_m4,EP,FT_rms_ind);
        ft_damage(ijk)=fqd;
    end    
%
end
%
pause(0.5);
progressbar(1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(fn,grms);
grid on;
xlabel('Natural Frequency (Hz)');
ylabel('Accel (GRMS)');
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
if(iu==1)
    out1=sprintf('Acceleration VRS  Q=%g  mass=%g lbm',Q,mass_orig);
else
    out1=sprintf('Acceleration VRS  Q=%g  mass=%g kg',Q,mass_orig);   
end
title(out1);
ymax=10^(ceil(log10(max(grms))));
ymin=10^(floor(log10(min(grms))));
ylim([ymin ymax]);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(fn,vrms);
grid on;
xlabel('Natural Frequency (Hz)');
if(iu==1)
    ylabel('Vel (in/sec RMS)');
    out1=sprintf('Velocity VRS  Q=%g  mass=%g lbm',Q,mass_orig);
else
    ylabel('Vel (m/sec RMS)');
    out1=sprintf('Velocity VRS  Q=%g  mass=%g kg',Q,mass_orig);    
end
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
title(out1);
ymax=10^(ceil(log10(max(vrms))));
ymin=10^(floor(log10(min(vrms))));
ylim([ymin ymax])
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
if(iu==2)
    drms=drms*1000;
end
plot(fn,drms);
grid on;
xlabel('Natural Frequency (Hz)');
if(iu==1)
    ylabel('Disp (in RMS)');
    out1=sprintf('Displacement VRS  Q=%g  mass=%g lbm',Q,mass_orig);
else
    ylabel('Disp (mm RMS)');
    out1=sprintf('Displacement VRS  Q=%g  mass=%g kg',Q,mass_orig);    
end
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
title(out1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(fn,FT_rms);
grid on;
xlabel('Natural Frequency (Hz)');
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
if(iu==1)
    ylabel('Transmitted Force (lbf RMS)');
    out1=sprintf('Transmitted Force VRS  Q=%g  mass=%g lbm',Q,mass_orig);
else
    ylabel('Transmitted Force (N RMS)');
    out1=sprintf('Transmitted Force VRS  Q=%g  mass=%g kg',Q,mass_orig);    
end
title(out1);
ymax=10^(ceil(log10(max(FT_rms))));
ymin=10^(floor(log10(min(FT_rms))));
ylim([ymin ymax])
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(ird==1)
    figure(fig_num);
    fig_num=fig_num+1;
    plot(fn,a_damage);
    grid on;
    xlabel('Natural Frequency (Hz)');
    ylabel('Relative Damage Index');
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
    if(iu==1)
        out1=sprintf('Acceleration Fatigue Damage Spectrum \n mass=%g lbm  Q=%g  b=%g',mass_orig,Q,bex);
    else
        out1=sprintf('Acceleration Fatigue Damage Spectrum \n mass=%g kg   Q=%g  b=%g',Q,mass_orig,Q,bex);        
    end
    title(out1); 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(fn,ft_damage);
    grid on;
    xlabel('Natural Frequency (Hz)');
    ylabel('Relative Damage Index');
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','on');
%
    if(iu==1)
        out1=sprintf('Transmitted Force Fatigue Damage Spectrum \n mass=%g lbm  Q=%g  b=%g',mass_orig,Q,bex);
    else
        out1=sprintf('Transmitted Force Fatigue Damage Spectrum \n mass=%g kg   Q=%g  b=%g',Q,mass_orig,Q,bex);        
    end
    title(out1);
%
    out10=sprintf('\n Fatigue damage is calculated from: (peak-valley)/2 \n');
    disp(out10);  
%
end