function L = gaspl(R,f,T,p,den)
%gaspl    Path loss due to atmosphere gaseous absorption
%   L = gaspl(R,F,T,P,DEN) returns the path loss, L (in dB), due to
%   atmosphere gaseous absorption. R is an length-M vector whose entries
%   represent the propagation distances (in meters). F is an length-N
%   vector whose entries represent the signal carrier frequency (in Hz). T
%   is a scalar indicating the temperature (in Celsius). P is a scalar
%   specifying the dry air pressure (in Pa) and DEN is a scalar specifying
%   the water vapour density (in g/m^3).
%
%   L is an MxN matrix whose elements represents the path loss of each
%   propagation path under corresponding frequency.
%
%   The ITU atmosphere gas model is valid between 1 to 1000 GHz.
%
%   % Examples:
%
%   % Example 1:
%   %   Compute the path loss of a 3 km link at 1 GHz due to atmosphere
%   %   gases. Assume the temperature is 15 degrees Celsius. The dry air
%   %   pressure is 1013 hPa and the water vapour density is 7.5 g/m^3.
%   
%   L = gaspl(3000,1e9,15,1013*1e2,7.5)
%
%   % Example 2:
%   %   Plot the attenuation model for atmosphere gases from 1 GHz to 1000
%   %   GHz. The dry air pressure is 1013 hPa and the water vapour density
%   %   is 7.5 g/m^3.
%
%   freq = (1:1000)*1e9;
%   L = gaspl(1e3,freq,15,101300,7.5);
%   semilogy(freq/1e9,L); grid on;
%   xlabel('Frequency (GHz)'); ylabel('Gas attenuation (dB/km)')
%
%   See also cranerainpl, fogpl, fspl, rainpl.

%   Copyright 2015-2016 The MathWorks, Inc.

% References
% [1] Recommendation ITU-R P.676-10, 2013
% [2] John Seybold, Introduction to RF Propagation, Wiley, 2005

%#codegen

validateattributes(R, {'double'}, {'nonnan','nonempty','finite','real', ...
    'nonnegative','vector'}, 'gaspl', 'R');
validateattributes(f, {'double'}, {'finite','nonempty','real', ...
    'vector','positive','>=',1e9,'<=',1000e9}, 'gaspl', 'F');
validateattributes(T, {'double'}, {'nonempty','finite','real', ...
    'scalar'}, 'gaspl', 'T');
validateattributes(p, {'double'}, {'nonnan','nonempty','real', ...
    'nonnegative','scalar','finite'}, 'gaspl', 'P');
validateattributes(den, {'double'}, {'finite','nonempty','real', ...
    'nonnegative','scalar'}, 'gaspl', 'DEN');

gamma = gasatt(f(:).',T,p,den);
Rkm = R(:)/1e3;
L = bsxfun(@times,Rkm,gamma);

end

function gamma = gasatt(f,Tc,p,den)

fGHz = f(:).'/1e9;
PhPa = p/100;
T = Tc + 273.15; % Convert C to K
theta = 300/T;
e = den*T/216.7;  % eq(4)

[Svec_oxy,Fmat_oxy] = dryaircoeff(fGHz,theta,PhPa,e);
[Svec_watvap,Fmat_watvap] = watvapcoeff(fGHz,theta,PhPa,e);

d = 5.6e-4*(PhPa+e)*theta^0.8;
NDdp = fGHz*PhPa*theta.^2.*(6.14e-5./(d.*(1+(fGHz./d).^2))+...
    1.4e-12*PhPa*theta^1.5./(1+1.9e-5.*fGHz.^1.5));

gammao = zeros(size(NDdp));
fidx = (fGHz<=118.750343);
gammao(fidx) = sum(bsxfun(@times,Svec_oxy,Fmat_oxy(:,fidx)));
fiidx = 38;  % pg 2
gammao(~fidx) = sum(bsxfun(@times,Svec_oxy(fiidx:end),Fmat_oxy(fiidx:end,~fidx)));
gammao = 0.1820*fGHz.*(gammao+NDdp);

gammaw = sum(bsxfun(@times,Svec_watvap,Fmat_watvap));
gammaw = 0.1820*fGHz.*gammaw;

gamma = gammao + gammaw;

end

function [Sivec,Fimat] = dryaircoeff(fGHz,theta,PhPa,e)

oxytab = [...
    % f0 ,a1 ,a2 ,a3 ,a4 ,a5 ,a6 
    50.474214 ,0.975 ,9.651 ,6.690 ,0.0 ,2.566 ,6.850 ;...
    50.987745 ,2.529 ,8.653 ,7.170 ,0.0 ,2.246 ,6.800 ;...
    51.503360 ,6.193 ,7.709 ,7.640 ,0.0 ,1.947 ,6.729 ;...
    52.021429 ,14.320 ,6.819 ,8.110 ,0.0 ,1.667 ,6.640 ;...
    52.542418 ,31.240 ,5.983 ,8.580 ,0.0 ,1.388 ,6.526 ;...
    53.066934 ,64.290 ,5.201 ,9.060 ,0.0 ,1.349 ,6.206 ;...
    53.595775 ,124.600 ,4.474 ,9.550 ,0.0 ,2.227 ,5.085 ;...
    54.130025 ,227.300 ,3.800 ,9.960 ,0.0 ,3.170 ,3.750 ;...
    54.671180 ,389.700 ,3.182 ,10.370 ,0.0 ,3.558 ,2.654 ;...
    55.221384 ,627.100 ,2.618 ,10.890 ,0.0 ,2.560 ,2.952 ;...
    55.783815 ,945.300 ,2.109 ,11.340 ,0.0 ,-1.172 ,6.135 ;...
    56.264774 ,543.400 ,0.014 ,17.030 ,0.0 ,3.525 ,-0.978 ;...
    56.363399 ,1331.800 ,1.654 ,11.890 ,0.0 ,-2.378 ,6.547 ;...
    56.968211 ,1746.600 ,1.255 ,12.230 ,0.0 ,-3.545 ,6.451 ;...
    57.612486 ,2120.100 ,0.910 ,12.620 ,0.0 ,-5.416 ,6.056 ;...
    58.323877 ,2363.700 ,0.621 ,12.950 ,0.0 ,-1.932 ,0.436 ;...
    58.446588 ,1442.100 ,0.083 ,14.910 ,0.0 ,6.768 ,-1.273 ;...
    59.164204 ,2379.900 ,0.387 ,13.530 ,0.0 ,-6.561 ,2.309 ;...
    59.590983 ,2090.700 ,0.207 ,14.080 ,0.0 ,6.957 ,-0.776 ;...
    60.306056 ,2103.400 ,0.207 ,14.150 ,0.0 ,-6.395 ,0.699 ;...
    60.434778 ,2438.000 ,0.386 ,13.390 ,0.0 ,6.342 ,-2.825 ;...
    61.150562 ,2479.500 ,0.621 ,12.920 ,0.0 ,1.014 ,-0.584 ;...
    61.800158 ,2275.900 ,0.910 ,12.630 ,0.0 ,5.014 ,-6.619 ;...
    62.411220 ,1915.400 ,1.255 ,12.170 ,0.0 ,3.029 ,-6.759 ;...
    62.486253 ,1503.000 ,0.083 ,15.130 ,0.0 ,-4.499 ,0.844 ;...
    62.997984 ,1490.200 ,1.654 ,11.740 ,0.0 ,1.856 ,-6.675 ;...
    63.568526 ,1078.000 ,2.108 ,11.340 ,0.0 ,0.658 ,-6.139 ;...
    64.127775 ,728.700 ,2.617 ,10.880 ,0.0 ,-3.036 ,-2.895 ;...
    64.678910 ,461.300 ,3.181 ,10.380 ,0.0 ,-3.968 ,-2.590 ;...
    65.224078 ,274.000 ,3.800 ,9.960 ,0.0 ,-3.528 ,-3.680 ;...
    65.764779 ,153.000 ,4.473 ,9.550 ,0.0 ,-2.548 ,-5.002 ;...
    66.302096 ,80.400 ,5.200 ,9.060 ,0.0 ,-1.660 ,-6.091 ;...
    66.836834 ,39.800 ,5.982 ,8.580 ,0.0 ,-1.680 ,-6.393 ;...
    67.369601 ,18.560 ,6.818 ,8.110 ,0.0 ,-1.956 ,-6.475 ;...
    67.900868 ,8.172 ,7.708 ,7.640 ,0.0 ,-2.216 ,-6.545 ;...
    68.431006 ,3.397 ,8.652 ,7.170 ,0.0 ,-2.492 ,-6.600 ;...
    68.960312 ,1.334 ,9.650 ,6.690 ,0.0 ,-2.773 ,-6.650 ;...
    118.750334 ,940.300 ,0.010 ,16.640 ,0.0 ,-0.439 ,0.079 ;...
    368.498246 ,67.400 ,0.048 ,16.400 ,0.0 ,0.000 ,0.000 ;...
    424.763020 ,637.700 ,0.044 ,16.400 ,0.0 ,0.000 ,0.000 ;...
    487.249273 ,237.400 ,0.049 ,16.000 ,0.0 ,0.000 ,0.000 ;...
    715.392902 ,98.100 ,0.145 ,16.000 ,0.0 ,0.000 ,0.000 ;...
    773.839490 ,572.300 ,0.141 ,16.200 ,0.0 ,0.000 ,0.000 ;...
    834.145546 ,183.100 ,0.145 ,14.700 ,0.0 ,0.000 ,0.000 ];

Sivec = oxytab(:,2)*1e-7*PhPa*theta^3.*exp(oxytab(:,3)*(1-theta)); % eq(3)
deltafvec = oxytab(:,4)*1e-4.*(PhPa*theta.^(0.8-oxytab(:,5))+1.1*e*theta); % eq(6a)
deltafvec = sqrt(deltafvec.^2+2.25e-6);  % eq(6b)
deltavec = (oxytab(:,6)+oxytab(:,7)*theta)*1e-4*(PhPa+e)*theta^0.8; % eq(7)

fdiff = bsxfun(@minus,oxytab(:,1),fGHz);
fsum = bsxfun(@plus,oxytab(:,1),fGHz);
Fimat = bsxfun(@minus,deltafvec,bsxfun(@times,deltavec,fdiff))./bsxfun(@plus,fdiff.^2,deltafvec.^2)+...
    bsxfun(@minus,deltafvec,bsxfun(@times,deltavec,fsum))./bsxfun(@plus,fsum.^2,deltafvec.^2);
Fimat = bsxfun(@rdivide,fGHz,oxytab(:,1)).*Fimat;  % eq(5)

end

function [Sivec,Fimat] = watvapcoeff(fGHz,theta,PhPa,e)

watvaptab = [...
    % f0 ,b1 ,b2 ,b3 ,b4 ,b5 ,b6 
    22.235080 ,0.1130 ,2.143 ,28.11 ,.69 ,4.800 ,1.00 ;...
    67.803960 ,0.0012 ,8.735 ,28.58 ,.69 ,4.930 ,.82 ;...
    119.995940 ,0.0008 ,8.356 ,29.48 ,.70 ,4.780 ,.79 ;...
    183.310091 ,2.4200 ,.668 ,30.50 ,.64 ,5.300 ,.85 ;...
    321.225644 ,0.0483 ,6.181 ,23.03 ,.67 ,4.690 ,.54 ;...
    325.152919 ,1.4990 ,1.540 ,27.83 ,.68 ,4.850 ,.74 ;...
    336.222601 ,0.0011 ,9.829 ,26.93 ,.69 ,4.740 ,.61 ;...
    380.197372 ,11.5200 ,1.048 ,28.73 ,.54 ,5.380 ,.89 ;...
    390.134508 ,0.0046 ,7.350 ,21.52 ,.63 ,4.810 ,.55 ;...
    437.346667 ,0.0650 ,5.050 ,18.45 ,.60 ,4.230 ,.48 ;...
    439.150812 ,0.9218 ,3.596 ,21.00 ,.63 ,4.290 ,.52 ;...
    443.018295 ,0.1976 ,5.050 ,18.60 ,.60 ,4.230 ,.50 ;...
    448.001075 ,10.3200 ,1.405 ,26.32 ,.66 ,4.840 ,.67 ;...
    470.888947 ,0.3297 ,3.599 ,21.52 ,.66 ,4.570 ,.65 ;...
    474.689127 ,1.2620 ,2.381 ,23.55 ,.65 ,4.650 ,.64 ;...
    488.491133 ,0.2520 ,2.853 ,26.02 ,.69 ,5.040 ,.72 ;...
    503.568532 ,0.0390 ,6.733 ,16.12 ,.61 ,3.980 ,.43 ;...
    504.482692 ,0.0130 ,6.733 ,16.12 ,.61 ,4.010 ,.45 ;...
    547.676440 ,9.7010 ,.114 ,26.00 ,.70 ,4.500 ,1.00 ;...
    552.020960 ,14.7700 ,.114 ,26.00 ,.70 ,4.500 ,1.00 ;...
    556.936002 ,487.4000 ,.159 ,32.10 ,.69 ,4.110 ,1.00 ;...
    620.700807 ,5.0120 ,2.200 ,24.38 ,.71 ,4.680 ,.68 ;...
    645.866155 ,0.0713 ,8.580 ,18.00 ,.60 ,4.000 ,.50 ;...
    658.005280 ,0.3022 ,7.820 ,32.10 ,.69 ,4.140 ,1.00 ;...
    752.033227 ,239.6000 ,.396 ,30.60 ,.68 ,4.090 ,.84 ;...
    841.053973 ,0.0140 ,8.180 ,15.90 ,.33 ,5.760 ,.45 ;...
    859.962313 ,0.1472 ,7.989 ,30.60 ,.68 ,4.090 ,.84 ;...
    899.306675 ,0.0605 ,7.917 ,29.85 ,.68 ,4.530 ,.90 ;...
    902.616173 ,0.0426 ,8.432 ,28.65 ,.70 ,5.100 ,.95 ;...
    906.207325 ,0.1876 ,5.111 ,24.08 ,.70 ,4.700 ,.53 ;...
    916.171582 ,8.3400 ,1.442 ,26.70 ,.70 ,4.780 ,.78 ;...
    923.118427 ,0.0869 ,10.220 ,29.00 ,.70 ,5.000 ,.80 ;...
    970.315022 ,8.9720 ,1.920 ,25.50 ,.64 ,4.940 ,.67 ;...
    987.926764 ,132.1000 ,.258 ,29.85 ,.68 ,4.550 ,.90 ;...
    1780.000000 ,22300.0000 ,.952 ,176.20 ,.50 ,30.500 ,5.00 ];

Sivec = watvaptab(:,2)*1e-1*e*theta^3.5.*exp(watvaptab(:,3)*(1-theta)); % eq(3)
deltafvec = watvaptab(:,4)*1e-4.*(PhPa*theta.^watvaptab(:,5)+watvaptab(:,6)*e.*theta.^watvaptab(:,7)); % eq(6a)
deltafvec = 0.535*deltafvec+sqrt(0.217*deltafvec.^2+2.1316e-12*watvaptab(:,1).^2/theta);  % eq(6b)
deltavec = zeros(size(deltafvec)); % eq(7)

fdiff = bsxfun(@minus,watvaptab(:,1),fGHz);
fsum = bsxfun(@plus,watvaptab(:,1),fGHz);
Fimat = bsxfun(@minus,deltafvec,bsxfun(@times,deltavec,fdiff))./bsxfun(@plus,fdiff.^2,deltafvec.^2)+...
    bsxfun(@minus,deltafvec,bsxfun(@times,deltavec,fsum))./bsxfun(@plus,fsum.^2,deltafvec.^2);
Fimat = bsxfun(@rdivide,fGHz,watvaptab(:,1)).*Fimat;  % eq(5)

end

% [EOF]