function IP=INNPOL(D,v1,v2,v3,borg,epsA,D0)
set(gcf,'Pointer','Watch');

[m,n]=size(D);
I=[];
mD = max(max(abs(D)));
for i=1:m
    if min(D(i,:))/max(abs(D(i,:)))>-0.1 && max(D(i,:))/mD>1.E-4
        I=[I, i];
    end
end
D=D(I,:);

[m,n]=size(D);
e=ones(n,1);
%project rows of D into plane spanned by v2 and v3
if borg==1
    n=0;
    for i=1:m
            n=n+1;
            B(n,1)=D(i,:)*v2';
            B(n,2)=D(i,:)*v3';
%         
    end
else
    n=0;
    for i=1:m

            n=n+1;
            B(i,1)=D(i,:)*v2'/(D(i,:)*v1');
            B(i,2)=D(i,:)*v3'/(D(i,:)*v1');
    end
end

%compute the convex hull of the points
K=convhull(B(:,1),B(:,2));
[m,n]=size(K);
IP=zeros(m-1,2);
for i=1:(m-1)
    IP(i,1:2)=B(K(i),1:2);
end

im=m-1;

IP=IP(1:im,:);
test=0;
i=1;
n=zeros(2,1);
while test<im
   n(1)=IP(mod(i-2,im)+1,2)-IP(mod(i,im)+1,2); 
   n(2)=IP(mod(i,im)+1,1)-IP(mod(i-2,im)+1,1);
   n=n/norm(n,2);
   c=IP(mod(i,im)+1,1:2)*n;
   d=IP(i,1:2)*n-c;
   if abs(d)<1.0e-14
       test=0;
       if i<im
            IP(i:im-1,:)=IP(i+1:im,:);
       end
       im=im-1;
       i=i-1;
       IP=IP(1:im,:);
    else
        v=IP(mod(i-2,im)+1,1:2)-IP(i,1:2);
        w=IP(i,1:2)-IP(mod(i,im)+1,1:2);
        phi=acos(v*w'/(norm(v,2)*norm(w,2)));
        if phi<1.0e-3
            IP(i:im-1,:)=IP(i+1:im,:);
            im=im-1;
            i=i-1;
            test=0;
        else
            test=test+1;
        end
    end
   i=mod(i,im)+1;
end
IP=IP(1:im,:);  

[x_max,x_max_ind]=max(IP(:,1));
k=find(IP(:,1)==x_max);
[numb_x_max,h]=size(k);
x_max_ind=k(1);
if numb_x_max>1
    for i=2:numb_x_max
        if IP(x_max_ind,2)>IP(k(i),2)
            x_max_ind=k(i);
        end
    end
end
if x_max_ind>1
    h=IP(1:(x_max_ind-1),1:2);
    IP(1:(im+1-x_max_ind),:)=IP(x_max_ind:im,:);
    IP((im-x_max_ind+2):im,:)=h;
end  
    
    
    
    
    
