clear all;
clc ;
close all;

T_0 = 293.15;
T_01 = 273.16 ;
T = 20 + 273.15; 
p_s0 = 1;
F = logspace(1,6);

hr = 0;
psat = p_s0*10^(-6.8346*(T_01/T)^1.261 + 4.6151);
h = p_s0*(hr)*(psat/p_s0);
F_rO = 1/p_s0*(24 + 4.04*10^4*h*(0.02+h)/(0.391+h));
F_rN = 1/p_s0*(T_0/T)^(1/2)*( 9 + 280*h*exp(-4.17*((T_0/T)^(1/3)-1)) );
alpha_ps_0 = 100*F.^2./p_s0.*( 1.84*10^(-11)*(T/T_0)^(1/2)...
    + (T/T_0)^(-5/2)*(0.01278*exp(-2239.1/T)./(F_rO + F.^2/F_rO)...
    + 0.1068*exp(-3352/T)./(F_rN + F.^2/F_rN) ) );
a_ps_0 = alpha_ps_0*20/log(10);

hr = 10; 
psat = p_s0*10^(-6.8346*(T_01/T)^1.261 + 4.6151);
h = p_s0*(hr)*(psat/p_s0);
F_rO = 1/p_s0*(24 + 4.04*10^4*h*(0.02+h)/(0.391+h));
F_rN = 1/p_s0*(T_0/T)^(1/2)*( 9 + 280*h*exp(-4.17*((T_0/T)^(1/3)-1)) );
alpha_ps_10 = 100*F.^2./p_s0.*( 1.84*10^(-11)*(T/T_0)^(1/2) ...
    + (T/T_0)^(-5/2)*(0.01278*exp(-2239.1/T)./(F_rO + F.^2/F_rO) ...
    + 0.1068*exp(-3352/T)./(F_rN + F.^2/F_rN) ) );
a_ps_10 = alpha_ps_10*20/log(10);

hr = 20; 
psat = p_s0*10^(-6.8346*(T_01/T)^1.261 + 4.6151);
h = p_s0*(hr)*(psat/p_s0);
F_rO = 1/p_s0*(24 + 4.04*10^4*h*(0.02+h)/(0.391+h));
F_rN = 1/p_s0*(T_0/T)^(1/2)*( 9 + 280*h*exp(-4.17*((T_0/T)^(1/3)-1)) );
alpha_ps_20 = 100*F.^2./p_s0.*( 1.84*10^(-11)*(T/T_0)^(1/2) ...
    + (T/T_0)^(-5/2)*(0.01278*exp(-2239.1/T)./(F_rO + F.^2/F_rO) ...
    + 0.1068*exp(-3352/T)./(F_rN + F.^2/F_rN) ) );
a_ps_20 = alpha_ps_20*20/log(10);

hr = 40; 
psat = p_s0*10^(-6.8346*(T_01/T)^1.261 + 4.6151);
h = p_s0*(hr)*(psat/p_s0);
F_rO = 1/p_s0*(24 + 4.04*10^4*h*(0.02+h)/(0.391+h));
F_rN = 1/p_s0*(T_0/T)^(1/2)*( 9 + 280*h*exp(-4.17*((T_0/T)^(1/3)-1)) );
alpha_ps_40 = 100*F.^2./p_s0.*( 1.84*10^(-11)*(T/T_0)^(1/2) ...
    + (T/T_0)^(-5/2)*(0.01278*exp(-2239.1/T)./(F_rO + F.^2/F_rO) ...
    + 0.1068*exp(-3352/T)./(F_rN + F.^2/F_rN) ) );
a_ps_40 = alpha_ps_40*20/log(10);

hr = 60; 
psat = p_s0*10^(-6.8346*(T_01/T)^1.261 + 4.6151);
h = p_s0*(hr)*(psat/p_s0);
F_rO = 1/p_s0*(24 + 4.04*10^4*h*(0.02+h)/(0.391+h));
F_rN = 1/p_s0*(T_0/T)^(1/2)*( 9 + 280*h*exp(-4.17*((T_0/T)^(1/3)-1)) );
alpha_ps_60 = 100*F.^2./p_s0.*( 1.84*10^(-11)*(T/T_0)^(1/2) ...
    + (T/T_0)^(-5/2)*(0.01278*exp(-2239.1/T)./(F_rO + F.^2/F_rO) ...
    + 0.1068*exp(-3352/T)./(F_rN + F.^2/F_rN) ) );
a_ps_60 = alpha_ps_60*20/log(10);

hr = 80; 
psat = p_s0*10^(-6.8346*(T_01/T)^1.261 + 4.6151);
h = p_s0*(hr)*(psat/p_s0);
F_rO = 1/p_s0*(24 + 4.04*10^4*h*(0.02+h)/(0.391+h));
F_rN = 1/p_s0*(T_0/T)^(1/2)*( 9 + 280*h*exp(-4.17*((T_0/T)^(1/3)-1)) );
alpha_ps_80 = 100*F.^2./p_s0.*( 1.84*10^(-11)*(T/T_0)^(1/2) ...
    + (T/T_0)^(-5/2)*(0.01278*exp(-2239.1/T)./(F_rO + F.^2/F_rO) ...
    + 0.1068*exp(-3352/T)./(F_rN + F.^2/F_rN) ) );
a_ps_80 = alpha_ps_80*20/log(10);

hr = 100;
psat = p_s0*10^(-6.8346*(T_01/T)^1.261 + 4.6151);
h = p_s0*(hr)*(psat/p_s0);
F_rO = 1/p_s0*(24 + 4.04*10^4*h*(0.02+h)/(0.391+h));
F_rN = 1/p_s0*(T_0/T)^(1/2)*( 9 + 280*h*exp(-4.17*((T_0/T)^(1/3)-1)) );
alpha_ps_100 = 100*F.^2./p_s0.*( 1.84*10^(-11)*(T/T_0)^(1/2) ...
    + (T/T_0)^(-5/2)*(0.01278*exp(-2239.1/T)./(F_rO + F.^2/F_rO)...
    + 0.1068*exp(-3352/T)./(F_rN + F.^2/F_rN) ) );
a_ps_100 = alpha_ps_100*20/log(10);

psvg = figure (1);
loglog(F,a_ps_0, F,a_ps_10, F,a_ps_20, F,a_ps_40,...
    F,a_ps_60, F,a_ps_80, F,a_ps_100,'LineWidth',1);
xlabel({'f / p_s [Hz/atm]';'Frequency/pressure'},'FontSize',10,...
    'FontWeight','normal','FontName','Times New Roman');
ylabel('Absorption coefficient/pressure a / p_s [dB/100 m atm]',...
    'FontSize',10,'FontWeight','normal','FontName','Times New Roman');
title({'Sound absorption coefficient per atmosphere for air at 20°C ';...
    'according to relative humidity per atmosphere'},...
    'FontSize',10,'FontWeight','bold','FontName','Times New Roman')
hleg = legend('   0','   10',...
    '   20','   40','   60',...
    '   80','   100');
v = get(hleg,'title');
set(v,'string',{'h_r / p_s [%/atm]'},'FontName','Times New Roman','FontSize',10,...
    'BackgroundColor', 'white','EdgeColor','white','HorizontalAlignment','center');
set(hleg,'Location','SouthEast','EdgeColor','black')
axis([1e1 1e6 1e-3 1e4]);
grid on;
set(gca,'gridlinestyle','-');
set(gca,'MinorGridLineStyle','-')
plot2svg('Absorption_coefficient.svg',psvg);