c23456789012345678901234567890123456789012345678901234567890123456789012
c        1         2         3         4         5         6         7

      subroutine sleep2
        implicit none
        integer i, j
        do i=1, 10000
          j = j + 1
        enddo
      end subroutine

      double precision function norma(x, y)
        implicit none
        double precision x, y
        norma = sqrt(x*x + y*y)
        return
      end function

c     (Ax, Ay) - (Bx, By) - pierwszy odcinek
c     (Cx, Cy) - (Dx, Dy) - drugi odcinek
c     (Ox, Oy) - punkt przecięcia
c     Zwracana wartość: .true. - nastąpiło przecięcie, .false. - nie
      logical function przeciecie_odcinkow(Ax, Ay, Bx, By, Cx, Cy,
     +Dx, Dy, Ox, Oy)
       implicit none
       double precision Ax, Ay, Bx, By, Cx, Cy, Dx, Dy, Ox, Oy
       double precision BAx, Bay, DCx, DCy, a, b, c, d, e, f
       double precision W, eps, x, y
       parameter (eps=1e-15)
       BAx = Bx - Ax
       BAy = By - Ay
       DCx = Dx - Cx
       DCy = Dy - Cy
c układ równań:
c   ax + by = c
c   dx + ey = f
       a = BAy
       b = -BAx
       c = DCy
       d = -DCx
       e = a*Ax + b*Ay
       f = d*Cx + e*Cy
       W = a*e - b*d
       if (W.lt.eps) then
         przeciecie_odcinkow = .false.
         return
       endif
      x = (c*e - b*f) / W
      y = (a*f - c*d) / W
c a = rzut (x, y) na wektor BA
c b = rzut (x, y) na wektor DC
c jeżeli (x, y) na odcinku BA, to 0 <= a <= 1, podobnie dla DC
      a = BAx*x + BAy*y
      b = DCx*x + DCy*y
      if (a.lt.0.or.a.gt.1.or.b.lt.0.or.b.gt.1) then
        przeciecie_odcinkow = .false.
        return
      endif
      przeciecie_odcinkow = .true.
      Ox = x
      Oy = y
      return
      end function

      subroutine zastosuj_brzegi(i)
        implicit none
        double precision p(0:1), p_nowe(0:1), v(0:1), O(0:1)
        double precision banda(0:2, 0:1), norm(0:2, 0:1), ee, tmp
        common p, p_nowe, v, banda, norm, ee
        integer i, i2
        logical trafiony, przeciecie_odcinkow

        if (i.lt.2) then
          i2 = i + 1
        else
          i2 = 0
        endif
        trafiony = przeciecie_odcinkow(banda(i, 0), banda(i, 1),
     +  banda(i2, 0), banda(i2, 1), p(0), p(1), p_nowe(0), p_nowe(1),
     +  O(0), O(1))
        if (trafiony) then
c zatrzymaj punkt na punkcie zderzenia
          p_nowe(0) = O(0)
          p_nowe(1) = O(1)
c zmień skłądową prostopadłą prędkości
          tmp = v(0)*norm(i, 0) + v(1)*norm(i, 1)
          v(0) = v(0) - ee*tmp*norm(i, 0)
          v(1) = v(1) - ee*tmp*norm(i, 1)
        endif
      end subroutine


      program bilard
        implicit none
        double precision p(0:1), p_nowe(0:1), v(0:1), banda(0:2, 0:1)
        double precision norm(0:2, 0:1), min_v, dt, t, e, tmp, ee
        double precision norma
        integer i
        parameter (min_v = 0.1)
        parameter (dt = 0.01)
        common p, p_nowe, v, banda, norm, ee

        e = 0
        ee = e + 1.0

        p(0) = 0.0
        p(1) = 0.0

        v(0) = 3.0
        v(1) = 0.0

        banda(0, 0) = 0.0
        banda(0, 1) = 3.0
        banda(1, 0) = -2.0
        banda(1, 1) = -1.0
        banda(2, 0) = 2.0
        banda(2, 1) = -1.0


        norm(0, 0) = banda(0, 1) - banda(1, 1)
        norm(0, 1) = banda(1, 0) - banda(0, 0)
        tmp = norma(norm(0, 0), norm(0, 1))
        norm(0, 0) = norm(0, 0) / tmp
        norm(0, 1) = norm(0, 1) / tmp
        norm(1, 0) = banda(1, 1) - banda(2, 1)
        norm(1, 1) = banda(2, 0) - banda(1, 0)
        tmp = norma(norm(1, 0), norm(1, 1))
        norm(1, 0) = norm(1, 0) / tmp
        norm(1, 1) = norm(1, 1) / tmp
        norm(2, 0) = banda(2, 1) - banda(0, 1)
        norm(2, 1) = banda(0, 0) - banda(2, 0)
        tmp = norma(norm(2, 0), norm(2, 1))
        norm(2, 0) = norm(2, 0) / tmp
        norm(2, 1) = norm(2, 1) / tmp

        t = 0.0
        print *, "czas polozenie_x polozenie_y predkosc"
        do
          p_nowe(0) = p(0) + v(0)*dt
          p_nowe(1) = p(1) + v(1)*dt
          do i=1,3
            call zastosuj_brzegi(i)
          enddo
          p(0) = p_nowe(0)
          p(1) = p_nowe(1)
          tmp = norma(v(0), v(1))
          print *, t, p(0), p(1), tmp
          if( min_v.gt.tmp ) exit
          t = t + dt
c          call sleep(1)
        enddo
      end program