function results = variogram_sim_plot(N,varargin)

%VARIOGRAM_PROFILE_PLOT simulates grid pattern of normal variates and plots the 
%                       variogram

%Written by: TONY E. SMITH, 1/14/01

%Functions called: VARIOGRAM_NEW

%
%INPUTS:   N = grid size
%    varargin = optional value of bandwidth 
%
%OUTPUTS: 

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

%OUTPUT (i) DAT = (bX3)-matrix with DAT(i)=(Di,Vi,Ni),i=1:b
%                Di = distance value for bin i
%                Vi = variogram value for bin i
%                Ni = number of point in bin i
%       (ii) data = NxN matrix of generated data

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 < 2)
   
   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 simulation

data = randn(N,N);

% Alter corner data

%data(1,1) = abs(data(1,1));
%data(1,N) = abs(data(1,N));
%data(N,1) = abs(data(N,1));
%data(N,N) = abs(data(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

V = DAT(:,2);

%OUTPUT

results.DAT = DAT;

results.data = data;

%Now plot output

plot(D,V,'.','MarkerSize',15);

hold on;

ONE = ones(length(D),1);

plot(D,ONE,'r-');