%function teta = dTeta(dist)
% dTeta - calculates coefficient of "good data"
%	teta = dTeta(dist)
%	dist - matrix of distances between points
%	teta - matrix of coeficients, same size as dist
%	

global dist teta

tic;
DiffTol = 0.05;
MinStd  = 5;
Diff2Std= 20;
StdTol  = 4;
AsymCoef = 0.1;
MeanDist = mean(mean(dist));
WindLen  = 50;
n = size(dist,2);
h = n;
k = size(dist,1)/n;
teta = ones(size(dist));
x1 = zeros(size(dist,1),1); x2 = zeros(size(dist,1),1);
q1 = zeros(size(dist,1)+WindLen-1,1); q2 = q1;
a = 1;

for i = 1:n-1, 
disp(i);
    for j = i+1:n,
	x1 = dSpikeEliminate(dExtract( dist, i,j));	s1 = std(x1);
	ddx = diff(diff(x1)); 	sddx1 = std(dSpikeEliminate(ddx));
	DdxLimit = min([ (s1*sddx1/10)  50 ]);
	if( s1 < MinStd  | sddx1 > Diff2Std ),
	    q1 = zeros(size(q1));
	else,
	    % Looking for discontinuities (second derivative high) and marking the 
	    % points immediately after as bad
	    q1 = 1./(1+a*conv( dTriac(ddx,DdxLimit)/DdxLimit, ...
		 	exp(-(0:WindLen-1)/(WindLen/5))) );
	    % Looking for points that are too far from mean and marking them 
	    q1 = q1(1:k) .* (abs(x1-mean(x1)) < 3*s1);
	end; %if

	x2 = dSpikeEliminate(dExtract( dist, j,i));	s2 = std(x2);
	ddx = diff(diff(x2)); 	sddx2 = std(dSpikeEliminate(ddx));
	DdxLimit = min([ (s2*sddx2/10)  50 ]);
	if( s2 < MinStd  | sddx2 > Diff2Std  ),
	    q2 = zeros(size(q2));
	else,
	    q2 = 1./(1+a*conv( dTriac(ddx,DdxLimit)/DdxLimit, ...
		 	exp(-(0:WindLen-1)/(WindLen/5))) );
	    q2 = q2(1:k) .* (abs(x2-mean(x2)) < 3*s2);
	end; %if
	if (s1 >= MinStd & s2 >= MinStd & sddx1<=Diff2Std & sddx2<=Diff2Std  ),
	    % x1 and x2 should be equal, finding points to far from each other
	    q = abs(x1-x2) ./ (x1+x2) < AsymCoef; 	
	    q1 = q1 .* (q & (q2>0.2));  q2 = q2 .* (q & (q1>0.2));
	end;
	dist((0:k-1)*h+i, j ) = x1;
	dist((0:k-1)*h+j, i ) = x2;
	teta((0:k-1)*h+i, j ) = q1(1:k);
	teta((0:k-1)*h+j, i ) = q2(1:k);

end; end;
toc

clear AsymCoef MeanDist MinStd WindLen DiffTol StdTol Diff2Std DdxLimit         
clear   ans  i j   q2 ddx  x1 x2 n s1 s2  q sddx1 sddx2 q1   

%save distteta dist teta h w k 

%j=j+1; subplot(211); plot(dExtract( dist, i,j)); subplot(212); plot(dExtract( teta, i,j));