function [h,T]=wvprmtrs(sp,fr,df)
% WVPRMTRS Computes wave parameters from Doppler spectra (WERA) 
%
%   [h,T]=wvprmtrs(sp,fr,df) computes the root mean square wave height (h)
%   in meters and the corresponding mean wave period (T) in seconds from 
%   M doppler spectra given in matrix sp. Input variables fr and df are the
%   radar (radio) frequency and the frequency interval of the doppler
%   spectrum ordinates respectively (both in Hz). Matrix sp contains M rows
%   corresponding to M Doppler spectra and 512 columns corresponding to the 
%   Doppler frequencies of each spectrum (standard WERA format). 
%
%   The significant wave height is obtained as Hs=4*h.
%
%   The algorithm is based on the work by Barrick, D. E., "Extraction of 
%   wave parameters from measured HF radar sea-echo Doppler spectra," Radio 
%   Science, 12-3, pp 415-424, 1977. 
%

% Last revision May 23, 2002. RRH. Parts of the code originally developed by 
% HCG. 

% Version 1, Revision 1. 


% DEFINE CONSTANTS 

g=9.81;							% gravity constant 
c=3e8;							% speed of electromagnetic waves 
kr=2*pi*fr/c;					% radar wave number 
[M,N]=size(sp);
fB=sqrt(g*fr/(pi*c));		    % Bragg frequency 
fBil=N/2-round(fB/df);	        % index of left (negative) 1st order (Bragg) peak
fBir=N/2+round(fB/df);	        % index of right (positive) 1st order (Bragg) peak

% FIND HIGHEST FIRST ORDER PEAKS 

[pk,np]=max(sp');				% all peaks 
p=find(np<123);					% correction for peaks below -1.0 Hz 
if ~isempty(p)
    [pk(p),np(p)]=max(sp(p,[fBil-40:fBil+40])');
	np(p)=np(p)+fBil-41;        % (peak within 0.30 Hz from theoretical value)
end
q=find(np>389);				    % correction for peaks above 1.0 Hz 
if ~isempty(q)
	[pk(q),np(q)]=max(sp(q,[fBir-40:fBir+40])');
	np(q)=np(q)+fBir-41;
end

% COMPUTE NOISE LEVEL (SIDE OF HIGHEST PEAK ONLY)

noise=NaN*ones(M,1);
inpl=find(np<N/2);		    	% peaks in left (negative) side 
noise(inpl)=mean(sp(inpl,1:27)');
                                % mean of left-most 0.20 Hz of the spectrum

inpr=find(np>=N/2);	    		% peaks in right (positive) side
noise(inpr)=mean(sp(inpr,485:512)');
                                % mean of right-most 0.20 Hz of the spectrum

% FIND CURRENT SHIFT OF FIRST ORDER PEAKS 

np_s(inpl)=np(inpl)-fBil;	% left side 
np_s(inpr)=np(inpr)-fBir;   % right side 

% CORRECTION OF SPECTRAL ORDINATES 

sp_s=zeros(M,N);
   
% Spectra whose highest peak is in the left (negative) side of doppler 
% frequencies 
   
for k=1:length(inpl)
   inplk=inpl(k);
   if np_s(inplk)>0
      sp_s(inplk,N+1-np_s(inplk):N)=ones(1,abs(np_s(inplk)))*noise(inplk);
      sp_s(inplk,1:N-np_s(inplk))=sp(inplk,1+np_s(inplk):N);
      						% shift ordinates np_s units to the left
	else
      sp_s(inplk,1:-np_s(inplk))=ones(1,abs(np_s(inplk)))*noise(inplk);
      sp_s(inplk,1-np_s(inplk):N)=sp(inplk,1:N+np_s(inplk));
      						% shift ordinates np_s units to the right
	end
end

% Spectra whose highest peak is in the right (positive) side of doppler 
% frequencies 
   
for k=1:length(inpr)
   inprk=inpr(k);
   if np_s(inprk)>0
      sp_s(inprk,1:N-np_s(inprk))=sp(inprk,1+np_s(inprk):N);
      						% shift ordinates np_s units to the left
   else
      sp_s(inprk,1:-np_s(inprk))=ones(1,abs(np_s(inprk)))*noise(inprk);
      sp_s(inprk,1-np_s(inprk):N)=sp(inprk,1:N+np_s(inprk));
      						% shift ordinates np_s units to the right
   end
end

% FIND LIMITS OF FIRST AND SECOND ORDER REGIONS 
   
% Spectra whose highest peak is in the left (negative) side of doppler 
% frequencies 

inpl_l=zeros(length(inpl),1);
inpl_r=zeros(length(inpl),1);

for k=1:length(inpl)
    inplk=inpl(k);
	difl=diff(log10(sp_s(inplk,fBil-21:fBil-1)));
	dummyl=find(difl(1:19)<0 & difl(2:20)>0)+1;
	[dummy,indxl]=min(sp_s(inplk,fBil-22+dummyl));
						    % find position of troughs in left side of 
                            % peak (within 20 frequency intervals)
    difr=diff(log10(sp_s(inplk,fBil+1:fBil+21)));
    dummyr=find(difr(1:19)<0 & difr(2:20)>0)+1;
    [dummy,indxr]=min(sp_s(inplk,fBil+dummyr));
						    % find position of troughs in right side of 
						    % peak (within 20 frequency intervals)
	 if isempty(dummyl) | isempty(dummyr)
                            % check for flat sides (partial or total)
% %         warning(['Unable to calculate first order peak limits in spectrum ',...
%             int2str(inpl(k))]);
      	inpl_l(inplk)=fBil-21;
							% set limit to minimum point in interval
	    inpl_r(inplk)=fBil+21;
						    % set limit to maximum point in interval
   else
     	inpl_l(inplk)=fBil-22+dummyl(indxl);
							% set limit at minimum trough in interval
        inpl_r(inplk)=fBil+dummyr(indxr);
							% set limit at minimum trough in interval
   end
   
end

% Spectra whose highest peak is in the right (positive) side of doppler 
% frequencies 

inpr_l=zeros(length(inpr),1);
inpr_r=zeros(length(inpr),1);

for k=1:length(inpr)
	inprk=inpr(k);
	difr=diff(log10(sp_s(inprk,fBir+1:fBir+21)));
    dummyr=find(difr(1:19)<0 & difr(2:20)>0)+1;
						    % find position of troughs in right side of 
                            % peak (within 20 frequency intervals)
    [dummy,indxr]=min(sp_s(inprk,fBir+dummyr));
   
    difl=diff(log10(sp_s(inprk,fBir-21:fBir-1)));
    dummyl=find(difl(1:19)<0 & difl(2:20)>0)+1;
						    % find position of troughs in left side of 
						    % peak (within 20 frequency intervals)
    [dummy,indxl]=min(sp_s(inprk,fBir-22+dummyl));

	if isempty(dummyl) | isempty(dummyr)
                            % check for flat sides (partial or total)
        warning(['Unable to calculate first order peak limits in spectrum ',...
           int2str(inpr(k))]);
	    inpr_r(inprk)=fBir+21;
							% set limit to minimum point in interval
    	inpr_l(inprk)=fBir-21;
							% set limit to maximum point in interval
   else
	    inpr_r(inprk)=fBir+dummyr(indxr);
							% set limit to minimum trough in interval
    	inpr_l(inprk)=fBir-22+dummyl(indxl);
							% set limit to minimum trough in interval
	end
   
end

% INTEGRATION OF RADAR CROSS SECTIONS (SIDE OF HIGHEST PEAK ONLY) 

wfn=weightf(fr,N,df);		% Barrick's weighting function 

Ircs1=zeros(M,1);
Ircs2=zeros(M,1);
Ircs2p=zeros(M,1);

% Spectra with peaks in left side 

for k=1:length(inpl)
   inplk=inpl(k);
   Ircs1(inplk)=sum(sp_s(inplk,inpl_l(inplk):inpl_r(inplk))-noise(inplk));
							% integral of 1st order radar cross section    
   rcs2=(sp_s(inplk,:)-noise(inplk))./wfn;
   rcs2(inpl_l(inplk):inpl_r(inplk))=0.0;
   Ircs2(inplk)=sum(rcs2(1:floor(N/2)));
   							% integral of weighted 2nd order rcs
   f1=((1:inpl_l(inplk))-ceil(N/2)).*df;
   Ircs2p1=sum(rcs2(1:inpl_l(inplk)))./...
      sum(abs(f1/(-fB)-1).*rcs2(1:inpl_l(inplk)));
   f2=((inpl_r(inplk):floor(N/2))-ceil(N/2)).*df;
   Ircs2p2=sum(rcs2(inpl_r(inplk):floor(N/2)))./...
      sum(abs(f2/(-fB)-1).*rcs2(inpl_r(inplk):floor(N/2)));
   Ircs2p(inplk)=(Ircs2p1+Ircs2p2)/2;
							% integral factor for period estimation 
end

% Spectra with peaks in right side 

for k=1:length(inpr)
   inprk=inpr(k);
   Ircs1(inprk)=sum(sp_s(inprk,inpr_l(inprk):inpr_r(inprk))-noise(inprk));
							% integral of 1st order radar cross section    
   rcs2=(sp_s(inprk,:)-noise(inprk))./wfn;
   rcs2(inpr_l(inprk):inpr_r(inprk))=0.0;
   Ircs2(inprk)=sum(rcs2(ceil(N/2)+1:N));
							% integral of weighted 2nd order rcs
   f1=((inpr_r(inprk):N)-ceil(N/2)).*df;
   Ircs2p1=sum(rcs2(inpr_r(inprk):N))./...
      sum(abs(f1/fB-1).*rcs2(inpr_r(inprk):N));
   f2=((ceil(N/2)+1:inpr_l(inprk))-ceil(N/2)).*df;
   Ircs2p2=sum(rcs2(ceil(N/2)+1:inpr_l(inprk)))./...
      sum(abs(f2/fB-1).*rcs2(ceil(N/2)+1:inpr_l(inprk)));
   Ircs2p(inprk)=(Ircs2p1+Ircs2p2)/2;
							% integral factor for period estimation 
end

% COMPUTATION OF WAVE PARAMETERS 

h=1/kr.*sqrt(2*Ircs2./Ircs1);
h=real(h);					% set high noise spectra estimates to zero
T=1/fB.*Ircs2p;