function h_plot = sed_normplot(x,epdf)

%     	function coef = sed_normplot(x,epdf)
% NORMPLOT Displays a normal probability plot.
%	This version to be used to make a probability plot of
%	DISCRETE data, where x is the data, and epdf is the estimated
%	probability density function of the data; ie it is the
%	fraction of data withing each x-bin.
%
%	changes made to normplot.m by ckh9c 1/14/94 
%
%	H_PLOT = SED_NORMPLOT(X,epdf) makes a normal probability plot of the  
%	data in X. For matrix, X, SED_NORMPLOT displays a plot for each column.
%	H_PLOT is a handle to the plotted lines.
%
%       epdf(i) should be the FRACTION (0-1) of sediment within size class
%       phi(i) to phi(i-1). 
%	
%	The purpose of a normal probability plot is to graphically assess
%	whether the data in X could come from a normal distribution. If the
%	data are normal the plot will be linear. Other distribution types 
%	will introduce curvature in the plot.  

%	Copyright (c) 1993 by The MathWorks, Inc.
%	$Revision: 1.3 $  $Date: 1993/10/01 14:21:23 $


[nx,mx]=size(x);
[ne, me] = size(epdf);
if ((ne~=nx) & (me~=mx))
  disp('error - reorient your input variables')
  disp('nx ~=ne or mx~=me')
  return
end
if (ne==nx)
%% reorient them for them.
   epdf = epdf'
   x    = x';
end
[n,m]=size(x);
[ne, me] = size(epdf);

[sx i]= sort(-x);
% want x in descending order.
sx=x(i); 

if me==1
  epdf=epdf(i);
else
  for jj=1:ne
    epdf(jj,:)=epdf(jj,i);
  end
end

minx  = min(sx(1,:));
maxx  = max(sx(n,:));
range = maxx-minx;

if range>0
%  minxaxis  = minx-0.025*range;
%  maxxaxis  = maxx+0.025*range;
   minxaxis  = floor(minx)-1;
   maxxaxis  = ceil(maxx);
else
  minxaxis  = minx - 1;
  maxxaxis  = maxx + 1;
end

%ecdf = [0.5./n:1./n:(n - 0.5)./n];
%	estimate pdf from the cdf. - -ckh
ecdf=cumsum(epdf')'
indNaN=find((ecdf< 0.0001)|(ecdf> 0.9999))

y  = norminv(ecdf,0,1);
y(indNaN)=NaN*ones(length(indNaN),1);

%minyaxis  = norminv(0.25 ./n,0,1);
%maxyaxis  = norminv((n-0.25) ./n,0,1);


%p     = [0.001 0.003 0.01 0.02 0.05 0.10 0.25 0.5...
%         0.75 0.90 0.95 0.98 0.99 0.997 0.999];
p     = [0.0001 0.001 0.003 0.01 0.02 0.05 0.10 0.16 0.25 0.5...
         0.75 0.84 0.90 0.95 0.98 0.99 0.997 0.999 0.9999];

%set y axis to entire plot.
minyaxis = norminv(min(p),0,1);
maxyaxis = norminv(max(p),0,1);


%label1= str2mat('0.001','0.003', '0.01','0.02','0.05','0.10','0.25','0.50');
%label2= str2mat('0.75','0.90','0.95','0.98','0.99','0.997', '0.999');
label1= str2mat('0.0001','0.001','0.003', '0.01','0.02','0.05','0.10',...
		'0.16','0.25','0.50');
label2= str2mat('0.75','0.84','0.90','0.95','0.98','0.99','0.997', '0.999', '0.9999');

label = [label1;label2];

tick  = norminv(p,0,1);

%if nargout > 2
%   h = plot(sx,y,'+',qx,qy,'-',mx,my,'-.',polyval(coef,my),my,'--');
%else 
    h_plot=plot (sx,y,'+-');
%   plot(sx,y,'+',qx,qy,'-',mx,my,'-.');
%   plot(sx,y,'+',qx,qy,'-',mx,my,'-.',...
%	polyval(coef,[minyaxis maxyaxis]),[minyaxis maxyaxis],'--');
%   text (median(x),-2.5,sprintf('Mean Phi= %4.2f',coef(2)))
%   text (median(x),-3,sprintf('Std Dev= %4.3f',coef(1)))
%end	

set(gca,'YTick',tick,'YTickLabels',label);
set(gca,'YLim',[minyaxis maxyaxis],'XLim',[minxaxis maxxaxis]);

set(gca,'XDir','rev')

set(gca,'xlim',[x(end)-1, x(1)+1]);
xlabel('\Phi size');
ylabel('Fraction Finer Than');
title('Normal Probability Plot For Sediment Distribution');

grid;