function speedd=hforce(level0,level0size,rmap,r);

global intmat

[M,N]=size(rmap);
speedd=zeros(level0size,1);
kernel=dh((-r:r)');
y=zeros(2*r+1,1);
z=zeros(4,4);
v2=[level0(2,1)-level0(1,1),level0(2,2)-level0(1,2)];
d2=sqrt(v2(1)^2+v2(2)^2);

for i=2:level0size-1,
	v1=v2;
	v2=[level0(i+1,1)-level0(i,1),level0(i+1,2)-level0(i,2)];
	d1=d2;
	d2=sqrt(v2(1)^2+v2(2)^2);
	v=v1*d2+v2*d1;
	slope=-v(1)/v(2);

	%midx=(level0(i+1,1)+level0(i-1,1))/2;
	%midy=(level0(i+1,2)+level0(i-1,2))/2;
	%if midx~=level0(i,1),
	%	slope=(midy-level0(i,2))/(midx-level0(i,1));

	%if (level0(i+1,2)~=level0(i-1,2)),
	%	slope=(level0(i-1,1)-level0(i+1,1))/...
	%		(level0(i+1,2)-level0(i-1,2));
		tmp=sqrt(1+slope^2);
		points=[(-r:r)/tmp+level0(i,1);(-r:r)/tmp*slope+level0(i,2)]';
    	%else
	%	points=[ones(1,2*r+1)*level0(i,1);(-r:r)+level0(i,2)]';
    	%end
	lownei=floor(points);
	t=points(:,1)-lownei(:,1);
	u=points(:,2)-lownei(:,2);
	for j=1:2*r+1
		z(:)=intmat(:,(lownei(j,2)-1)*N+lownei(j,1));
		y(j)=[1 t(j) t(j)^2 t(j)^3]*z*[1 u(j) u(j)^2 u(j)^3]';
	end
	%y=(1-t).*(1-u).*mgrad((lownei(:,2)-1)*N+lownei(:,1))+...
	%	t.*(1-u).*mgrad((lownei(:,2)-1)*N+lownei(:,1)+1)+...
	%	(1-t).*u.*mgrad(lownei(:,2)*N+lownei(:,1))+...
	%	t.*u.*mgrad(lownei(:,2)*N+lownei(:,1)+1);
	if rmap(lownei(1,1),lownei(1,2))==1,
		speedd(i)=sum(kernel.*y);
	else
		speedd(i)=-sum(kernel.*y);
	end
end
speedd(1)=speedd(2);
speedd(level0size)=speedd(level0size-1);


