function OUT = variogram_profile_plot(N,sim,varargin)

%VARIOGRAM_PROFILE_PLOT simulates many grid patterns using data_sim and plots the 
%                       mean and standard deviation of values at each distance

%Written by: TONY E. SMITH, 1/14/01
%
%INPUTS:   N = grid size
%        sim = number of simulated grids
%    varagin = optional value of bandwidth 
%
%OUTPUTS: 

%ON SCREEN: Plot of the means for each dist plus variance (in red)

%OUT = (bX3)-matrix with OUT(i)=(Di,Mi,Si),i=1:b
%                Di = distance value for bin i
%                Mi = sample mean variogram value at distance Di
%                Si = sample standard deviation at distance Di

n = N^2 ; %number of locations

%Make locations

data = zeros(N,N);

[X,Y,junk] = find(data == 0);


%Compute matrix, D, of pairwise distances

U = ones(n,1) ; %column vector

XX = X*U' - U*X' ;

YY = Y*U' - U*Y' ;

D = sqrt(XX.^2 + YY.^2) ;

D = tril(D) ;  %take lower triangular part


%Now reduce to nonzero distances and determine range


[SS(:,1),SS(:,2),SS(:,3)] = find(D) ;

clear('D');

	%The matrix SS has rows (i,j,d) containing the row,
   %column, and values of all positive cells of D 


Dmin = min(SS(:,3)) ;

Dmax = max(SS(:,3)) ;

%Set BandWidth if not already specified

if (nargin < 3)
   
   if (Dmin < Dmax/2)
   
   d = Dmax/2 ; %maximum relevant distance
   
	else
   
      error('Dmin >= Dmax/2') ; 
      
   end
   
else
         
   d = varargin{:};
   
   if (Dmax < d)  %check to see if d is too large
      
      d = Dmax ; 
      
   end
   
end


%Do first simulation

data = randn(N,N);

val = reshape(data,N^2,1);

M = [X,Y,val];

DAT = variogram_new(M,d);

D = DAT(:,1); %vector of distances for each bin

Sum = DAT(:,2);

SumSqr = DAT(:,2).^2;

%Now do the rest of the simulations

for i = 2:sim
   
   data = randn(N,N);

	val = reshape(data,N^2,1);

	M = [X,Y,val];

DAT = variogram_new(M,d);

Sum = Sum + DAT(:,2);

SumSqr = SumSqr + DAT(:,2).^2;

completed = 100*(i/sim)
   
end

%Now calculate output

A = Sum/sim;

S = sqrt((SumSqr/sim) - A.^2);

OUT = [D,A,S];

plot(D,A,'.');

hold on;

plot(D,S,'r.');