program turtle
	implicit NONE
	real PI
	parameter (PI=3.141593)
	real x_min /-3.0/, x_max /3.0/
	real y_min /-3.0/, y_max /3.0/
	real radius /0.02/	
	real step /0.02/
	real x1(2) / -1.0, 0.0/
	real x2(2) /  1.0, 0.0/
	real step_rotate /30.0/
	integer lunout /3/
	integer err /-1/
c function:
	logical reach
c local:
	integer status
	real x(2),new_x(2)
	real field(2)
	real theta
1	continue
c	call file_open('output', lunout, status)
	if (status.eq.err) goto 1	
c	write(*,*) 'set size 10 by 10'
c	write(*,*) 'set limit x ',x_min,x_max
c	write(*,*) 'set limit y ',y_min,y_max
	theta = 15.0
	do while (theta.lt.355.)
	  x(1) = x1(1) + radius* cos(theta*pi/180.)
	  x(2) = x1(2) + radius* sin(theta*pi/180.)
	  write(*,*) '# theta:',theta
	do while (.not.reach(x,x2,step))
	  call calc_field(x1,x2,x,field)
	  call move (x,step,field,new_x)
	  x(1) = new_x(1)
	  x(2) = new_x(2)
	  if ((x(1).gt.x_min) .and. (x(1).lt.x_max)) then
             if ((x(2).gt.y_min) .and. (x(2).lt.y_max)) then
	      write(*,*) x
	  endif
	endif
	enddo
	write (*,*) '# plot'
	theta = theta + step_rotate
	end do
	end
	subroutine calc_field(x1,x2,x,field)
	implicit NONE
	real x1(2)
	real x2(2)
	real x(2)
	real field(2)
	real r1,r2
	  r1=sqrt((x(1)-x1(1))**2 + (x(2)-x1(2))**2)
	  r2=sqrt((x(1)-x2(1))**2 + (x(2)-x2(2))**2)
	  field(1)=(x(1)-x1(1))/(r1**3) + (-1.)*(x(1)-x2(1))/(r2**3)
	  field(2)=(x(2)-x1(2))/(r1**3) + (-1.)*(x(2)-x2(2))/(r2**3)
	return
	end
	subroutine move(x,step,field,new_x)
	implicit NONE
	real x(2)
	real step
	real field(2)
	real new_x(2)
	real eps
	parameter (eps=1.0e-20)
	real norm
		norm = sqrt(field(1)**2+field(2)**2)
		if (norm.lt.eps) then
		  write (*,*) 'too small'
		  new_x(1) = x(1)
		  new_x(2) = x(2)
		else
		  new_x(1) = x(1)+step*field(1)/norm
		  new_x(2) = x(2)+step*field(2)/norm
		endif
	return
	end
	logical function reach(x,point,distance)
	implicit NONE
	real x(2)	
	real point(2)
	real distance
	real r
	  r = sqrt((x(1)-point(1))**2+(x(2)-point(2))**2)
	  if (r.lt.distance) then
	    reach = .true.
	  else
	    reach = .false.
	  endif
	end