%
% Example on epage 489, 
% 
% Problem on epage 517
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');

randn('seed',0); 
rand('seed',0); 

load bodmin
% Loads data in x and y and boundary in bodpoly.
% Get the Ghat function first and plot.
X = [x,y];
n = length(x);

% The G function is the nearest neighbor 
% distances for each event.
% Find the distances for all points.
dist = pdist(X);
% Convert to a matrix and put large 
% numbers on the diagonal.
D = diag(realmax*ones(1,n)) + squareform(dist);
% Find the smallest distances in each row or col.
mind = min(D);

% create a sampling of nearest neighbor distances: 
dw = (max(mind)-min(mind))/30; 
w  = min(mind):dw:max(mind);
nw = length(w);
ghat = zeros(1,nw);

% Now get the values for ghat.
for i = 1:nw
  ind = find(mind<=w(i));
  ghat(i) = length(ind);
end
ghat = ghat/n;


% compute the expression 12.12 from using our estimate of ghat:
% 
% (a poor approximation of our bounding area is given by) 
r = ( max(bodpoly(:,1))-min(bodpoly(:,1)) )*( max(bodpoly(:,2))-min(bodpoly(:,2)) );

lambda_hat = length(x)/r; 

isqt = -log( ( 1 - ghat + 1e-6 )/(lambda_hat*pi) ); 
isqt(find(isqt<0))=0;
C = sqrt( isqt ); 
figure; plot( w, C, 'x-' ); 
xlabel('Event-Event Distances - w'); ylabel('C(w)'); grid on; 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_8', 'epsc' ); 

clear functions;