function [y,xn]=interpm(xa,ya,x)
%function [y,(x)]=interpm(xa,ya, x or n)
%      a no-loop interpolation routine
%      x is the evenly sampler array. Optionally n is number of points

if length(x)==1
    n=x;
%                      generate x array
    dx=(xa(length(xa))-xa(1))/(n-1);
    x=xa(1):dx:xa(length(xa));
    xn=x;
else
    n=length(x);
    dx=(x(n)-x(1))/(n-1);
end

%           get xa and ya in ascending order in x
%[xa,i]=sort(xa);
%ya=ya(i);
%                    get shape right
[i,j]=size(x);
[i1,j1]=size(xa);
if (i-j)*(i1-j1)<0
         x=x';
end

ind=ceil((xa-x(1))/dx)+1;
a=find(ind<1);
if length(a)>0 
         ind(a)=[];
end
ibelow=zeros(1:length(x));
ibelow(ind)=1+length(a):length(xa);
ind=find(ibelow~=0);
ibelow(ind)=[ibelow(ind(1)) diff(ibelow(ind))];
ibelow=cumsum(ibelow(1:n));
a=find(ibelow<1);
if length(a)>0
         ibelow(a)=ones(a/a(1));
end
%..........................................   now upper side
ind=floor((xa-x(1))/dx)+1;
a=find(ind<1);
if length(a)>0 
         ind(a)=[];
end
iabove=zeros(1:length(x));
%ind=ind(length(ind):-1:1);
iabove(ind)=1+length(a):length(xa);
ind=find(iabove~=0);
iabove(ind)=[iabove(ind(1)) diff(iabove(ind))];
a=length(x):-1:1;
iabove(a)=cumsum(iabove(a));
iabove=iabove(1)+1-iabove(1:length(x));
a=find(iabove>length(xa));
if length(a)>0
         iabove(a)=ones(a/a(1))*length(xa);
end
y=ya(iabove);
a=find(iabove-ibelow~=0);
if length(a)>0
fac=(x(a)-xa(ibelow(a)))./(xa(iabove(a))-xa(ibelow(a)));
y(a)=ya(ibelow(a)).*(1-fac)+ya(iabove(a)).*fac;
end
