function [gamma, r, theta, diff2, numpair]=variogram2dL(x, dat, h, angles, r, theta, diff2)

%
%[gamma, r, theta, diff2]=variogram2d(x, dat, h, angles, r, theta, diff2)
%
%Computes the 2 dimensional (semi)variogram of the data in dat.  For a variogram in 2D,
%the lag distances are actually sectors of a disc defined in polar coordinates.  The
%sectors are specifed by h (radii) and angles (in radians...from -90 to 90 degrees).
%x is a column vector of positions i.e. x(i,:) is the position of the i'th data point and
%x should be n-by-d where n is the number of points and d=2 is the dimensionality of the data set.
%dat should be n-by-s where n is the number of points and s is the number of different parameters
%measured at each point.  The variogram gamma will then be a 3D array of size 
%length(h) by length(angles) by s and will contain the semivariance of each variable 
%at each sector.  r, theta, and diff2 are optional arguments used to
%calculate the variogram.  Passing them will greatly speed the computation.


%x:	column vector of positions
%dat:	data on which to calculate the variogram
%h:		vector of distances to bin h=100:100:5000;
%angles:	in radians from -90 to 90 degrees
%r:			
%theta:		
%diff2:		difference in value between pairs of every point


tic


[n,d]=size(x);
l=length(h);
l2=length(angles);
[tmp, s]=size(dat);
if(tmp ~= n)
	error(['x and dat must have the same number of rows']);
end;
if(d~=2)
	error('x must be 2-dimensional');
end;

numpairs=.5*n*(n-1); 	%number of pairs of points (i, j)--just an arithmetic series minus n since
					%we don't care about the difference between a point and itself

if nargin==4
	r=zeros(numpairs,1);
	theta=r;
	diff2=zeros(numpairs, s);
end
if nargin==5
	theta=zeros(numpairs,1);
	diff2=zeros(numpairs, s);
end;
if nargin==6
	if (length(r)~=numpairs | length(theta)~=numpairs)
		error('Incorrect number of pairs in r and/or theta');
	end
	diff2=zeros(numpairs, s);
end
if nargin==7
	[t1, t2]=size(diff2);
	if(t1~=numpairs & t2~=s)
		error('Incorrect number of pairs in diff2');
	end
end

tmp=1;
tmp2=1;
if nargin<7
	for k=1:n-1
		l1=n-k;
		l2=k+1;
		for i=1:l1;
			j=k+i;				%Louise added this based on coding in variogram.m
			if nargin==4
				r(tmp)=norm(x(i,:)-x(j,:));	%distance between pairs of points
			end;
			if nargin<=5
				if(r(tmp)~=0)
					theta(tmp)=asin((x(i,2)-x(j,2))/r(tmp));
				end;
			end;
			tmp=tmp+1;				%Louise added this
		end;
		I=find(r(tmp2:tmp2+l1-1)==0);
		diff2(tmp2:tmp2+l1-1,:)=(dat(1:l1,:)-dat(l2:n,:)).^2;
		diff2(tmp2-1+I,:)=ones(length(I), s)*NaN;
		tmp2=tmp2+l1;
	end;
end

l=length(h);
l2=length(angles);
if(h(1)==0)
	display('Variance at radius h=0 not computed');
else
	h=[0,h];
	l=l+1;
end;
gamma=NaN*ones(l-1, l2, s);
numpair=zeros(size(gamma));
if(min(diff(h))<=0)
	error(['Radii in h should be unique and increasing']);
end;

if(min(diff(angles))<=0)
	error(['angles should be unique and increasing']);
end;

for i=2:l;								
	I=find(r>h(i-1) & r<=h(i));
	if(length(I) > 0)
		for j=2:l2;
			%J=find(theta>angles(j-1) & theta<angles(j));  Louise replaced this line with the following line.
			J=find(theta>angles(j-1) & theta<=angles(j));
			if (length(J) > 0)
				Int=intersect(I, J);
				ll = length(Int);
				if (ll > 0)
					numpair(i-1,j,:)=ll; %Louise added this to count number of pairs used to calc gamma
					gamma(i-1,j,:)=sum(diff2(Int,:))/(2*ll);
				end;
			end;		
		end;
	end;
end;
bb = find(gamma == NaN);
display([int2str(length(bb)), ' bins out of ', int2str((l-1)*l), ' contained no data']);

toc


