c subroutine dwr(r,n,y,dr)
c r is r(n,3)
c  dr is the length of dwr part  
c  y(i) is between r(i-dr),r(i+dr) totally 2*dr+1 beads 2*dr segments
c the same as  ~/works/dwr/dwr.m

	subroutine dwr(r,n,y,dr)
        implicit double precision(a-h,o-z)
	parameter(maxa=500)
	double precision r(maxa,3),y(maxa),A(maxa,3)
	double precision ddwr,RR(3),A1(3),B1(3)
	integer n,dr,u,i,i2,i3,j,dr2
	
	do i=1,n
	y(i)=0.0
	enddo

c A(i) is the vector between r(i+1) and r(i)
	do i=1,n
 	 i2=i+1
	 if (i.eq.n) then 
	  i2=1
	 endif
	 do j=1,3
	  A(i,j)=r(i2,j)-r(i,j)
	 enddo
 	enddo
 
c u is the separation between two segments A(i) and A(i+u)
	 dr2=2*dr-1
	do u=2,dr2
	 do i=1,n
	  i2=i+u
	  if (i2.gt.n) then
	   i2=i2-n
	  endif
	  do j=1,3
	   RR(j)=r(i2,j)-r(i,j)
	   A1(j)=A(i,j)
	   B1(j)=A(i2,j)
  	  enddo
      	  ddwr=dwr_ab(A1,B1,RR)
c	  write(*,*) 'in dwr',A1(1),B1(1),RR(1)
c	  write(*,*) ddwr,i,i2
c	  pause
	   if (abs(ddwr).gt.3.0) then
		  ddwr=0.0
	   endif
	   do i3=i2+1-dr,i+dr
		   ii=i3
		   if (i3.le.0) then
	 	    ii=i3+n
		   endif
		   if (i3.gt.n) then
		    ii=i3-n
		   endif
	   y(ii)=y(ii)+ddwr
  	   enddo
	enddo
	enddo  

	return
	end