function h=lcplot(d,plotflag,ma,sa);
%LCPLOT Plots a level crossing spectrum (lc) 
%
% CALL:  lcplot(lc,plotflag,ma,sa)
%
%      lc  = a two column matrix with levels and number of upcrossings.
%
%plotflag = 1 plot level crossing  spectrum and
%           2 overplots it with the theoretical one for a Gaussian process
%             (default) 
%
%   ma,sa = mean and standard deviation of the process (default estimated from lc)
%
% See also: dat2lc, mm2lc, wnormpdf, polyfit

% Tested on: Matlab 5.3
% History:
% revised by jr 06.07.00. 
% - Division by 2 from 05.07.00 removed (line 42). 
% - sa -> sa^2 (line 66)
% revised by jr 05.07.00. Line 39: division by 2
% revised by pab 11.08.99
% removed the plot procedure from the old mm2cross function  
% added the option of overplotting of the theoretical one for a Gaussian process

if nargin<2|isempty(plotflag)
 plotflag=2
end
[cmax icmax]=max(d(:,2));% maximum crossing
if (nargin <4|isempty(sa))&plotflag==2, 
  % if stdev of x unknown then estimate it
  cros=d;
  indz=(cros(:,2)==0);
  lncros=cros; 
  logcros(indz,2)=inf; % this is done in order to avoid a warning message
  logcros(~indz,2)=-log(cros(~indz,2));
  logcmin =logcros(icmax,2);
  logcros(:,2)=sqrt(2*abs(logcros(:,2)-logcmin)); %sqrt(2*logcros)
  logcros(1:icmax,2)=2*logcros(icmax,2)-logcros(1:icmax,2); 
  p = polyfit(d(10:end-9,1),logcros(10:end-9,2),1); %least square fit
  
  sa=1/p(1); % estimated standard deviation of x
  
  %plot(d(:,1),logcros(:,2)),hold on
  %plot(d(:,1),polyval(p,d(:,1)),'r'),hold off, pause
end

if (nargin<3|isempty(ma))&plotflag==2, 
  ma=0; % the mean is apriori zero.
  %mu=p(2);%, d(icmax,1) % new estimated mean
  %mu=d(icmax,1); % maximum cr intensity should be at the mean
  % but this not always the case
end


%  clf
hh=stairs(d(:,1),d(:,2));
axis([-1.05*max(abs(d(:,1))) 1.05*max(abs(d(:,1))) -1 1.05*max(d(:,2))])
title('Number of upcrossings')
xlabel('level u')
%spclbox


% bar(d(:,1),d(:,2))
if plotflag==2
  y = wnormpdf(d(:,1),ma,sa^2);  
  y = cmax*y/y(icmax);%Normalization needed to overplot the crossingspectrum.
  np = get(gca,'NextPlot');    
  set(gca,'NextPlot','add')    
  hh1 = plot(d(:,1),y,'r--','LineWidth',1);% Plots density line over histogram.
  set(gca,'NextPlot',np) 
end

if nargout == 1 & plotflag==2
  h = [hh; hh1];
elseif nargout==1
  h=hh;
end