function [wfn,f,fB]=weightf(fr,N,df)
%WEIGHTF Calculates Barrick's weighting function for WERA (smooth)
%
%   wfn=weightf(fr,N,df) calculates the weighting function (Barrick's) to 
%   account for the contribution of wave-wave interactions in the second order 
%   peaks of doppler spectra. The frequency of the radar (in Hz) is given by 
%   fr. N is the number of doppler frequencies to be calculated (even) and df 
%   is the frequency interval of the components. 
%
%   [wfn,f]=weightf(fr,N,df) same as before, but giving also the frequency 
%   range as output. 
%
%   [wfn,f,fB]=weightf(fr,N,df) same as option above, but giving also the
%   Bragg frequency as output.
%

% Original version by SFH, 1996. Update and debugging by RRH, 2000. 
% Code optimization by RRH, April 10, 2002. 

% Definition of constant parameters 

g=9.81;						% gravity constant 
c=3e8;						% speed of electromagnetic waves
fB=sqrt(g*fr/(pi*c));	    % Bragg frequency
f=(1:(N/2))*df;	        	% frequency range 

% Definition of Barrick's weighting function. 

x1=[0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.66 0.7 0.8 0.9 1];
x1=fB*[x1 1.1 1.15 1.2 1.3 1.4 sqrt(2)];
y1=[1e4 8e2 65 18 6.5 3.3 1.8 1.7 1.74 1.98 2.6 2.8];
y1=[y1 2.65 2.60 2.65 2.9 4.7 5.3];
x2=fB*[sqrt(2) 1.49 1.5 1.55 1.6 2^(3/4)];
y2=[5.3 1.98 1.99 3 7.4 190];
x3=fB*[2^(3/4) 1.69 1.7 1.75 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5];
y3=[190 16 9 5.8 5.4 5.9 6.8 8.5 11 14.8 20 27];
pp4=polyfit(x3(11:12),log(y3(11:12)),1);
							% fit for frequency values > 1.23 Hz

% Separate frequency range of into 4 categories (corresponding to smooth 
% sections of Barrick's weighting function). 

indx1=find(f>=0 & f<=sqrt(2)*fB);
indx2=find(f<=2^(3/4)*fB & f>sqrt(2)*fB);
indx3=find(f<=2.4*fB & f>2^(3/4)*fB);
indx4=find(f>2.4*fB);

f1=f(indx1);
f2=f(indx2);
f3=f(indx3);
f4=f(indx4);

% Find ordinates of frequency values within each category (cubic 
% interpolation). 

wf1=exp(spline(x1,log(y1),abs(f1)));
wf2=exp(spline(x2,log(y2),abs(f2)));
wf3=exp(spline(x3,log(y3),abs(f3)));
wf4=exp(polyval(pp4,abs(f4)));

wfn=[wf1 wf2 wf3 wf4];
indx=length(wfn):-1:1;
wfn=[wfn(indx),wfn];

if nargout>1
   f=[-f(indx),f];
end