program aaa
use dflib
implicit none
integer i, j, k, m, counterBall, counter, mask(1000), balls, check, steps, next, cyrcle, linewidth, speedcolor, color(1000), resultspeed(250), graphsize
real x(1000), y(1000), way(1000), position, lx, ly, lx2, ly2, eee, d, d1, d2, wx1, wx, wy1, wy, tmp, speed(1000), eps, minspeed, angle, maxspeed, Macsvell, tmpspeed(1000), pi, addx, addy
type (xycoord) xy
call setvieworg(400,100,xy)

pi = 3.14159
balls = 1000 ! kolichestvo sharov
linewidth = 10 ! tolshina linii
cyrcle = 180 ! scorost' vrasheniya
next = 5 ! scorost' poyavleniya
eps = 0.0001 ! tochnost' integrala
minspeed = 0
maxspeed = 3
angle = pi / 12 ! ugol v radianax
graphsize = 250

counterBall = 1

1 continue

!call sleepqq(10)

if ((int(counter/next) + 1).ge.counterBall .and. counterBall.le.balls) then
    call random_seed()
    call random_number(tmp)
    way(counterBall) = tmp * pi * 32 !32 casuse more random.
    mask(counterBall) = 1
    speed(counterBall) = Macsvell(eps, minspeed, maxspeed)
    color(counterBall) = speedcolor(speed(counterBall), minspeed, maxspeed)
    counterBall = counterBall + 1
endif

call clearscreen($gclearscreen)
eee = -counter
eee = eee / cyrcle
wx = int(150*sin(eee+angle))
wy = int(150*cos(eee+angle))
wx1 = int(150*sin(eee-angle))
wy1 = int(150*cos(eee-angle))
d=setcolorrgb(#ff0000)
d2=ellipse ($GFILLINTERIOR, -linewidth,-linewidth,500+linewidth,500+linewidth)
d=setcolor(0)
d2=ellipse ($GFILLINTERIOR, 0,0,500,500)
d1=setcolor(12)
d1=pie ($GFILLINTERIOR, 100-linewidth,100-linewidth,400+linewidth,400+linewidth,250+wy1,250+wx1,250+wy,250+wx)
d1=setcolor (0)
d1=ellipse($GFILLINTERIOR, 100,100,400,400)
d=setcolor(9)
d=ellipse($GFILLINTERIOR, 245,245,255,255)

check = 0

do i = 1, balls

    if (mask(i).eq.1) then
        check = check + 1
        tmpspeed(check) = speed(i)
        eee = counter
        eee = eee / cyrcle
        steps = counter - i * next + next
        x(i) = cos(way(i)) * speed(i) * steps
        y(i) = sin(way(i)) * speed(i) * steps
        lx = x(i) + 250
        ly = y(i) + 250
        k = getpixel(lx, ly)
        !write(*, *) (x(i) ** 2 + y(i) ** 2), (250 ** 2)
        !call sleepqq(100)
        if (k.eq.12 .or. ((x(i) ** 2 + y(i) ** 2).ge.(250 ** 2))) then
            if (lx.lt.(400+linewidth) .and. lx.gt.(100-linewidth) .and. ly.gt.(100-linewidth) .and. ly.lt.(400+linewidth)) then
                mask(i) = 0
            else
                mask(i) = 2
                lx2 = cos(way(i)) * steps * speed(i)
                ly2 = sin(way(i)) * steps * speed(i)
                x(i) = lx2*cos(-eee) + ly2*sin(-eee)
                y(i) = ly2*cos(-eee) - lx2*sin(-eee)
            endif
        else
            d = setcolorrgb(color(i))
            d = ellipse($GFILLINTERIOR, lx-5, ly-5, lx+5, ly+5)
        endif
    endif

    if (mask(i).eq.2) then
        eee = counter
        eee = eee / cyrcle
        lx2 = x(i)
        ly2 = y(i)
        lx = lx2*cos(eee) + ly2*sin(eee) + 250
        ly = ly2*cos(eee) - lx2*sin(eee) + 250
        d = setcolor(10)
        d = ellipse($GFILLINTERIOR, lx-5, ly-1, lx+5, ly+1)
    endif
enddo

if (check.eq.0 .and. counterBall.ge.balls) then
    goto 2
endif

write(*,*) (tmpspeed(i),i=1,check)

counter = counter + 1
goto 1

2 continue

check = 0

do i = 1, balls
    if (mask(i).eq.2) then
        check = check + 1
    endif
    j = int((speed(i) - minspeed) * graphsize / (maxspeed - minspeed))
    resultspeed(j) = resultspeed(j) + 1
enddo

m = 0
do i = 1, graphsize
    m = max(m, resultspeed(i))
enddo

k = int(graphsize / (maxspeed - minspeed))
k = 16
do i = 1, graphsize
    do j = 0, resultspeed(i) * k
        d = setcolorrgb(Z"ffff" - Z"100" * int(255 * j / (m * k)))
        d = setpixel(-graphsize + i, 500 - j)
        d = setpixel(-graphsize + i + 1, 500 - j)
    enddo
enddo

write(*, *) check, balls

end



real function F(x)
implicit none
real x

F = 4*3.14159*(3.14159)**(1.5)*exp(-x**2)*x**2

return
end

real function integ(eps, l, r)
implicit none
real l, r, F, eps
integer num, i

integ = 0
num = (r-l) / eps
do i = 1, num
    integ = integ + eps * F(l + i * eps)
enddo

return
end

real function searchinteg(eps, l, r, value)
implicit none
real l, r, F, eps, integ, value
integer num, i

integ = 0
num = (r-l) / eps
do i = 1, num
    integ = integ + eps * F(l + i * eps)
    if (integ.ge.value) then
        searchinteg = l + i * eps - eps / 2
        return
    endif
enddo
searchinteg = r
return
end

real function Macsvell(eps, mins, maxs)
implicit none
integer m
real rand, eps, max, integ, mid, tmp, mins, maxs, searchinteg

call random_seed()
call random_number(rand)
max = integ(eps, mins, maxs)
rand = rand * max
Macsvell = searchinteg(eps, mins, maxs, rand)

return
end

integer function speedcolor(speed, min, max)
real speed, min, max

speedcolor = Z"440000" + Z"10001" * int(speed*160 / (max - min))

return
end