function [arcs,rhos]=curvature(xs)
%computes curvature function rhos
%for coordinate matrix xs

l=length(xs(:,1));
%discretize perimeter into arcs of approximately length dl
dl=25;arcs=[];

%extend perimeter by one period
xs=[xs;xs];t=1;

	while t<=l
   lengths=zeros(1,l);
   	for j=1:l
      lengths(j)=norm(xs(t+j,:)-xs(t,:));
   	end
   
      if (max(lengths)>dl)    
         [z,k]=max(lengths>=dl);
    		xlo=xs(t+k-1,:)-xs(t,:);xhi=xs(t+k,:)-xs(t,:);
   		l1=dl-lengths(k-1);l2=lengths(k)-dl;
   		c1=l2/(l1+l2);c2=l1/(l1+l2);
   		arc=c1*(xs(t+k,:)-xs(t,:))  +c2*(xs(t+k-1,:)-xs(t,:));
         arcs=[arcs;arc];
         t=t+k;
        else 
      end
   end
   
   %now get discretized curvatures rhos
   narcs=length(arcs(:,1));
   rhos=zeros(1,narcs-1);
   arcs=[arcs;arcs(1,:)];
   for t=1:narcs-1
      a=arcs(t,:);a=a./norm(a);
      b=arcs(t+1,:);b=b./norm(b);
      x=a*b';
      y=a(1)*b(2)-a(2)*b(1);
      r=sqrt(x^2+y^2);
      rhos(t)=(x>=0)*asin(y/r)+(x<0)*(pi-asin(y/r));
      
     % rhos(t)=2*rhos(t)/(norm(arcs(t,:))+norm(arcs(t+1,:)));
      
   end