function [box zero range eigvec]=make_box4(points)
points=min(points,20*ones(size(points,1),size(points,2)));
fi=8;
a0=0;
%OK
N=size(points,1);
dim=size(points,2);
num_of_points=round(N^(1/dim));

meanval=mean(points,1);
eps=0.1;
[eigvec,eigval]=eig((points-ones(N,1)*meanval)'*(points-ones(N,1)*meanval));
range=max(2*sqrt(diag(eigval)/N),eps);

%plot(points(:,1), points(:,2),'r.' ); hold on

pom=zeros(4,2);
p=1;
%figure
while p
    for i=1:4
        pom(i,:)=(range.*[(-1)^i; (-1)^floor((i-1)/2)])'*eigvec'+meanval;
    end
    %plot(pom(:,1),pom(:,2), 'r+'); hold on
    if min(pom(:))<0
        [a b]= min(pom(:,2));
        [a c]= max(pom(:,2));
        a1=(pom(b,1)-pom(c,1))>0;
        range=range/1.05;
        if range>eps
            if a1~=a0
                fi=2*fi;
            end
            a0=a1;
            eigvec=[cos(pi/fi) (-1)^a1*sin(pi/fi);(-1)^(a1+1)*sin(pi/fi) cos(pi/fi)]*eigvec;
        end
    else 
        p=0;
    end
end


%for i=1:2
 %   pom=eigvec(:,i)*range(i)+meanval';
  %  plot([meanval(1) pom(1)] , [meanval(2) pom(2)]);
%end


box=zeros(N,dim);
pointer=-range;
step=2*range/(num_of_points-1);
jump=[1 num_of_points]; 
for i=1:N
    for k=1:dim
        box(i,k)=pointer(k)+step(k)*mod(floor((i-1)/jump(k)),num_of_points);
    end
end  

box=(box)*eigvec'+ones(N,1)*meanval;
range=ones(dim,1)./(2*range);
zero=box(1,:);
end