program hs
  implicit none
  integer, parameter          :: n     = 10000
  double precision, parameter :: rad   = 0.25d0/sqrt(dble(n))
  integer, parameter          :: nstop = 100000
  integer            :: i, j, nrejet ! iterateurs
  logical            :: accept       ! condition d'acceptance du tirage
  double precision   :: x(n), y(n)   ! les mailles du reseau
  double precision   :: d2, a, b     ! variables intermediaires pour le tirage
  double precision   :: d2min        ! distance min entre 2 points au carre

  ! pas besoin de prendre des racines inutiles à chaque tirage et chaque test !
  d2min = (2.d0*rad)**2

  ! Le premier point... car il en faut bien un ;-) !
  call random_number(a)
  call random_number(b)
  x(1)=a
  y(1)=b
  print *, x(1), y(1), rad, 1, 0

  ! Les autres points
  do i=2,n
    nrejet = 0               ! On initialise un compteur quelconque ...
                             !    ... pour ne pas boucler indéfiniement
lm: do                       ! On commence un boucle pour tirer un point valide
      call random_number(a)  ! On tire x ...
      call random_number(b)  !           et y ...
      accept = .true.        !                et on suppose que c'est ok
lt:   do j=1,i-1             ! ET ON TESTE SUR TOUS LES POINTS DEJA VALIDES
        d2 = (x(j)-a)**2 + (y(j)-b)**2
        if ( d2 < d2min ) then
          accept = .false.   ! si on a un et un seul overlap, c'est mauvais !
          exit lt            !    => pas la peine de tester les autres points
        endif
      enddo lt
      if ( accept ) exit lm  ! si on a accept qui est encore true ici : YES !
                             !          ... et on peut sortir de la boucle lm
      nrejet = nrejet+1      ! sinon, on a un rejet supplémentaire et on boucle...
                             ! pas trop quand meme...
      if ( nrejet > nstop ) stop 'trop de rejets successifs...' 
    enddo lm
    ! si on en est arrive la sans STOP, on a notre point 'i' 
    x(i) = a
    y(i) = b
    print *, x(i), y(i), rad, i, nrejet
  enddo
!  do i=1,n
!    print *, x(i), y(i), rad, i
!  enddo

end program hs