Corona virus simulator

General FreeBASIC programming questions.
albert
Posts: 5916
Joined: Sep 28, 2006 2:41
Location: California, USA

Re: Corona virus simulator

Postby albert » Mar 26, 2020 2:19

They need to take white cells from people that have recovered the fastest , and implant them into new infected cases...

Make a human serum like a horse serum...
integer
Posts: 391
Joined: Feb 01, 2007 16:54
Location: usa

Re: Corona virus simulator

Postby integer » Mar 26, 2020 3:22

@dodicat
THUMBS UP!
UEZ
Posts: 627
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Postby UEZ » Mar 26, 2020 7:30

badidea wrote:
UEZ wrote:Alive means not dead.^^ The pixel is infected (red) but still alive and with a probability of 5% the pixel will die. Alive = infected (not dead) + cured / inoculated.

But this comment is confusing:
Dim As Ubyte status, zone 'status: alive = 0, infected = 1, cured = 2, dead = 3, inoculated = 4
On infection you switch form 0 to 1. So 0 must mean 'initial' or something I think.

To my 'forked' version, I added a 'day counter' and a 'infected by gunman (patient 0)' counter.
With patient 0 at the same speed as the others, he infects 10 to 25 others (depending on his actual speed and start position).
[code]dim shared as integer iGunMan, infectedbyGunMan = 0


Another speed optimisation is to set an end-of-sickness-date, and compare against that. So that you don't have to do Particles(i).days += 0.0075 for all infected every loop. This part could even be done 1 once every 10 loops. Currently, in simulated time, it correspond to a heath check every ~10 minutes.


Yes, "alive" is not properly chosen - healthy is a better wording for that.

Your settings looks more realistic and it is a good idea to display the passed days. I will modify my last version accordingly. :-)

I assume that there are a lot more optimization possibilities in the code. One of my first ideas was to create zones and check infection there only to speed it up but I don't know whether it was a good idea...

dodicat wrote:Not so much a simulation but a message.

Very nice idea. Keep well and fit.
ShawnLG
Posts: 137
Joined: Dec 25, 2008 20:21

Re: Corona virus simulator

Postby ShawnLG » Mar 26, 2020 14:28

UEZ wrote:
badidea wrote:
UEZ wrote:I assume that there are a lot more optimization possibilities in the code. One of my first ideas was to create zones and check infection there only to speed it up but I don't know whether it was a good idea...


There is optimizations available such as partitioning the space into a tree datastructure.
https://en.wikipedia.org/wiki/Nearest_neighbor_search It would be a good idea if you wanted to simulate millions of persons.
grindstone
Posts: 752
Joined: May 05, 2015 5:35
Location: Germany

Re: Corona virus simulator

Postby grindstone » Mar 26, 2020 16:17

@dodicat: Nice work!
badidea
Posts: 2149
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Mar 26, 2020 21:39

UEZ wrote:One of my first ideas was to create zones and check infection there only to speed it up but I don't know whether it was a good idea...

Yes, you can do a 'pre-selection'. You know that a particle or person in de top-right corner can never interact with bottom-left. Spilt the map in a grid. Only check within a grid-section, and neighbouring grid-sections. Make it flexible, so you can test what grid size is optimal. I did this before, I'll look up the code...

I found the code I wrote 8 years ago. A particle collision simulation. The variable bruteForce determines whether to check all particle against each other (1), or with pre-sectioning in a 8 x 8 grid (0). Will take me some time to understand my own code. I should have put more comments in. There is a nice feature: You can zoom in and out with the mouse-wheel and drag the view with right-mouse-button pressed.'

Code: Select all

type timer_type
   dim as integer active
   dim as double total
   dim as double started
   dim as string name
   declare sub start()
   declare sub stop()
end type

sub timer_type.start()
   if (active = 0) then
      active = 1
      started = timer()
   end if
end sub

sub timer_type.stop()
   if (active = 1) then
      total += (timer() - started)
      active = 0
   end if
end sub

#define nTimers 3

dim as timer_type tt(nTimers-1)
tt(0).name = "0"
tt(1).name = "1"
tt(2).name = "2"

'------------------- Default types and numbers ------------------

type sgl2d
   as single x, y
end type

type int2d
   as integer x, y
end type

const as single pi = 3.141592654
'const as integer true = 0
'const as integer false = -1

'---------------------- Mouse stuff ----------------------

#DEFINE MOUSE_IDLE 0
#DEFINE MOUSE_POS_CHANGED 1
#DEFINE MOUSE_LB_PRESSED 2
#DEFINE MOUSE_LB_RELEASED 3
#DEFINE MOUSE_RB_PRESSED 4
#DEFINE MOUSE_RB_RELEASED 5
#DEFINE MOUSE_MB_PRESSED 6
#DEFINE MOUSE_MB_RELEASED 7
#DEFINE MOUSE_WHEEL_UP 8
#DEFINE MOUSE_WHEEL_DOWN 9

type mouseType
 pos as int2d
 posChange as int2d
 wheel as integer
 buttons as integer
 lb as integer 'left button
 rb as integer 'right button
 mb as integer 'middle button
end type

function handleMouse(byref mouse as mouseType) as integer
  static previous as mouseType
  dim as integer change = MOUSE_IDLE
  getmouse mouse.pos.x, mouse.pos.y, mouse.wheel, mouse.buttons
  if (mouse.buttons = -1) then
    mouse.lb = 0
    mouse.rb = 0
    mouse.mb = 0
    mouse.posChange.x = 0
    mouse.posChange.y = 0
  else
    mouse.lb = (mouse.buttons and 1)
    mouse.rb = (mouse.buttons shr 1) and 1
    mouse.mb = (mouse.buttons shr 2) and 1
    if (previous.pos.x <> mouse.pos.x or previous.pos.y <> mouse.pos.y) then
      change = MOUSE_POS_CHANGED
    end if
    mouse.posChange.x = mouse.pos.x - previous.pos.x
    mouse.posChange.y = mouse.pos.y - previous.pos.y
    if (previous.buttons <> mouse.buttons) then
      if (previous.lb = 0 and mouse.lb = 1) then change = MOUSE_LB_PRESSED
      if (previous.lb = 1 and mouse.lb = 0) then change = MOUSE_LB_RELEASED
      if (previous.rb = 0 and mouse.rb = 1) then change = MOUSE_RB_PRESSED
      if (previous.rb = 1 and mouse.rb = 0) then change = MOUSE_RB_RELEASED
      if (previous.mb = 0 and mouse.mb = 1) then change = MOUSE_MB_PRESSED
      if (previous.mb = 1 and mouse.mb = 0) then change = MOUSE_MB_RELEASED
    end if
    if (mouse.wheel > previous.wheel) then change = MOUSE_WHEEl_UP
    if (mouse.wheel < previous.wheel) then change = MOUSE_WHEEl_DOWN
    previous = mouse
  end if
  return change
end function

'---------------------- Scaled graphics ----------------------

type scaled_graphics_type
   as single scale' = 1 / pixel_size 'pixels / meter
   as int2d offset' = (scrn_w \ 2, scrn_h \ 2) 'offset in pixels
   as integer scrn_w, scrn_h
   declare sub init(sc as single, xOffset as integer, yOffset as integer)
   declare sub pset_xy(x as single, y as single, c as integer)
   declare sub pset(p as sgl2d, c as integer)
   declare sub circle_xy(x as single, y as single, r as single, c as integer)
   declare sub circle(p as sgl2d, r as single, c as integer)
   declare sub line_xy(x1 as single, y1 as single, x2 as single, y2 as single, c as integer)
   declare sub line(p1 as sgl2d, p2 as sgl2d, c as integer)
end type

sub scaled_graphics_type.init(sc as single, xoffs as integer, yoffs as integer)
   scale = sc
   offset.x = xoffs
   offset.y = yoffs
   screeninfo scrn_w, scrn_h
end sub

sub scaled_graphics_type.pset_xy(x as single, y as single, c as integer)
   dim as integer xInt = int(offset.x + x * scale)
   dim as integer yInt = int(offset.y + y * scale)
   .pset(xInt, scrn_h - yInt), c
end sub

sub scaled_graphics_type.pset(p as sgl2d, c as integer)
   dim as integer xInt = int(offset.x + p.x * scale)
   dim as integer yInt = int(offset.y + p.y * scale)
   .pset(xInt, scrn_h - yInt), c
end sub

sub scaled_graphics_type.circle_xy(x as single, y as single, r as single, c as integer)
   dim as integer xInt = int(offset.x + x * scale)
   dim as integer yInt = int(offset.y + y * scale)
   .circle(xInt, scrn_h - yInt), r * scale, c
end sub

sub scaled_graphics_type.circle(p as sgl2d, r as single, c as integer)
   dim as integer xInt = int(offset.x + p.x * scale)
   dim as integer yInt = int(offset.y + p.y * scale)
   .circle(xInt, scrn_h - yInt), r * scale, c
end sub

sub scaled_graphics_type.line_xy(x1 as single, y1 as single, x2 as single, y2 as single, c as integer)
   dim as integer xInt1 = int(offset.x + x1 * scale)
   dim as integer yInt1 = int(offset.y + y1 * scale)
   dim as integer xInt2 = int(offset.x + x2 * scale)
   dim as integer yInt2 = int(offset.y + y2 * scale)
   .line(xInt1, scrn_h - yInt1)-(xInt2, scrn_h - yInt2), c
end sub

sub scaled_graphics_type.line(p1 as sgl2d, p2 as sgl2d, c as integer)
   dim as integer xInt1 = int(offset.x + p1.x * scale)
   dim as integer yInt1 = int(offset.y + p1.y * scale)
   dim as integer xInt2 = int(offset.x + p2.x * scale)
   dim as integer yInt2 = int(offset.y + p2.y * scale)
   .line(xInt1, scrn_h - yInt1)-(xInt2, scrn_h - yInt2), c
end sub

'---------------------- Vector stuff ----------------------

function vector_add(v1 as sgl2d, v2 as sgl2d) as sgl2d
   dim as sgl2d v
   v.x = v1.x + v2.x
   v.y = v1.y + v2.y
   return v
end function

function vector_mul(v1 as sgl2d, mul as single) as sgl2d
   dim as sgl2d v
   v.x = v1.x * mul
   v.y = v1.y * mul
   return v
end function

'----------------------------------------------------------
'---------------------- Main program ----------------------
'----------------------------------------------------------

'screen stuff
const as single pixel_size = 6e-11 'm
const as integer scrn_w = 800
const as integer scrn_h = 600

'physical constants
const as single g = 9.81 'm/s^2
const atomicMass as double = 1.66e-27 'kg
const univGasConst as double = 8.314 'J/mol K
const amuArgon as double = 39.95 'no unit
const mArgonMol as double = amuArgon / 1000.0 'kg/mol

'stuff to play with
const as double kBall = 10 'N/m
const as integer nBalls = 900
const as integer nSquares = 2
const as integer nSectx = 8
const as integer nSecty = 8
const as integer maxCollisions = 200
const as integer bruteForce = 0

dim shared as scaled_graphics_type sg

'---------------------- Squares ----------------------

type line_type
   A as sgl2d
   dxBA as single
   dyBA as single
   lSquared as single
end type

type quad_type
   line(4-1) as line_type
   faceingOut as integer
end type

sub drawQuad(square as quad_type, c as integer)
   dim as integer i, j
   dim as sgl2d B
   for i = 0 to 3
      j = (i + 1) and 3
      B.x = square.line(i).A.x + square.line(i).dxBA
      B.y = square.line(i).A.y + square.line(i).dyBA
      sg.line(square.line(i).A, B, c)
      'sg.line(square.line(i).A, square.line(j).A, c)
   next
'   sg.line(square.line(0).A, square.line(1).A, c)
'   sg.line(square.line(1).A, square.line(2).A, c)
'   sg.line(square.line(2).A, square.line(3).A, c)
'   sg.line(square.line(3).A, square.line(0).A, c)
end sub

sub updateQuad(quad as quad_type)
   dim as integer i, j
   with quad
      for i = 0 to 3
         j = (i + 1) and 3
         .line(i).dxBA = .line(j).A.x - .line(i).A.x
         .line(i).dyBA = .line(j).A.y - .line(i).A.y
         .line(i).lSquared = .line(i).dxBA * .line(i).dxBA + .line(i).dyBA * .line(i).dyBA
      next
   end with
end sub

sub makeSquare(quad as quad_type, size as single)
   with quad
      .faceingOut = true 'still unused
      .line(0).A.x = -size / 2:   .line(0).A.y = -size / 2
      .line(1).A.x =  size / 2:   .line(1).A.y = -size / 2
      .line(2).A.x =  size / 2: .line(2).A.y =  size / 2
      .line(3).A.x = -size / 2:   .line(3).A.y =  size / 2
   end with
   updateQuad(quad)
end sub

sub makeDiamond(quad as quad_type, size as single)
   with quad
      .faceingOut = true 'still unused
      .line(0).A.x = 0:     .line(0).A.y = -size
      .line(1).A.x = size:  .line(1).A.y = 0
      .line(2).A.x = 0:     .line(2).A.y = size
      .line(3).A.x = -size: .line(3).A.y = 0
   end with
   updateQuad(quad)
end sub

'Returns false if not on line segment
function closestPointOnLineSegment(C as sgl2d, A as sgl2d, B as sgl2d, P as sgl2d) as integer
   'point C, line AB
   dim as single dxBA = B.x - A.x
   dim as single dyBA = B.y - A.y
   dim as single dxCA = C.x - A.x
   dim as single dyCA = C.y - A.y
   dim as single lSquared_AB = dxBA * dxBA + dyBA * dyBA '(Bx-Ax)^2 + (By-Ay)^2
   dim as single dotProduct_AC_AB = dxCA * dxBA + dyCA * dyBA '(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
   dim as single rFactor = dotProduct_AC_AB / lSquared_AB
   'dim as single determinant_BA_AC = -dyCA * dxBA - -dxCA * dyBA  '(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
   'dim as single sFactor = determinant_BA_AC / lSquared_AB
   'Test if projection of C is on line segment
   'r=0      P = A
   'r=1      P = B
   'r<0      P is on the backward extension of AB
   'r>1      P is on the forward extension of AB
   '0<r<1    P is interior to AB   
   if (rFactor < 0) then
      return false
   else
      if (rFactor > 1) then
         return false
      else
         P.x = A.x + rFactor * dxBA
         P.y = A.y + rFactor * dyBA
         return true
      end if
   end if
end function

'Returns false if not on line segment - OPTIMIZED
function cpols(C as sgl2d, lijn as line_type, byref P as sgl2d) as integer
   dim as single dxCA = C.x - lijn.A.x
   dim as single dyCA = C.y - lijn.A.y
   dim as single dotProduct_AC_AB = dxCA * lijn.dxBA + dyCA * lijn.dyBA '(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
   dim as single rFactor = dotProduct_AC_AB / lijn.lSquared
   if (rFactor < 0) then
      return false
   else
      if (rFactor > 1) then
         return false
      else
         P.x = lijn.A.x + rFactor * lijn.dxBA
         P.y = lijn.A.y + rFactor * lijn.dyBA
         return true
      end if
   end if
end function

'---------------------- Balls ----------------------

type ball_type
   p as sgl2d 'position [m]
   v as sgl2d 'velocity [m/s]
   a as sgl2d 'accelleration [m/s^2]
   F as sgl2d 'force [N]
   m as single 'mass [kg]
   r as single 'radius [m]
end type

'---------------------- Sections ----------------------

type section_type
   as integer nItems
   as integer index(nBalls-1)
   declare sub addBall(iBall as integer)
end type

sub section_type.addBall(iBall as integer)
   index(nItems) = iBall
   nItems += 1
end sub

'---------------------- Collision lists ----------------------

type colPairList_type
   as integer nItems
   as integer index1(maxCollisions-1)
   as integer index2(maxCollisions-1)
   declare function addPairIfNew(iBall1 as integer, iBall2 as integer) as integer
   declare sub printList()
end type

sub colPairList_type.printList()
   dim as integer i
   for i = 0 to nItems-1
      print index1(i), index2(i)
   next
end sub

function colPairList_type.addPairIfNew(iBall1 as integer, iBall2 as integer) as integer
   dim as integer i
   'check new pair
   for i = 0 to nItems-1
      if (iBall1 = index1(i)) and (iBall2 = index2(i)) then return false
   next
   'add pair
   if (nItems < maxCollisions) then
      index1(nItems) = iBall1
      index2(nItems) = iBall2
      nItems += 1
   else
      print "Error: colPairList_type.addPair: maxCollisions"
   end if
   return true
end function

'---------------------- Variables ----------------------
   
dim as ball_type ball(nBalls-1)
dim as quad_type square(nSquares-1)
dim shared as section_type section(nSectx-1, nSecty-1) 'shared because of size
dim as colPairList_type colPairList
dim as single t, dt 's
dim as integer i, j, k, i2, j2, nList, result, collisionCount
dim as integer xi, yi, xi1, xi2, yi1, yi2
dim as single xBall, yBall, rBall
dim as ball_type ptr pBall
dim as single dy, dx, dCenter, dEdge, F
dim as single Ekin, Eelas, Etot
dim as sgl2d P, A, B
dim as integer drawCounter = 0, loopCounter = 0
dim as single temperature = 300 'K
dim as single vAverage, alfa
dim as single square1size = 3e-8, xSectDiv, ySectDiv

dim as mousetype mouse
dim as integer mouseEvent
dim as integer mouseDrag = 0

dim as double startTime = timer

'---------------------- Main loop ----------------------

screenres scrn_w, scrn_h, 32
sg.init(1 / pixel_size, scrn_w \ 2, scrn_h \ 2)

makeSquare(square(0), square1size)
makeDiamond(square(1), square1size / 3)

vAverage = sqr((3 * univGasConst * temperature) / mArgonMol) ' m/s

'randomize(timer)
randomize(12345)
for i = 0 to nBalls-1
   with ball(i)
      .m = amuArgon * atomicMass 'kg
      .r = 188 * 1e-12 '188 pm = Van der Waals radius
      .p.x = ((i mod 30) / 15 - 1.0) * 1.2e-8
      .p.y = ((i \ 30) / 15 - 1.0) * 1.2e-8
      alfa = rnd(1) * 2 * pi
      .v.x = cos(alfa) * vAverage
      .v.y = sin(alfa) * vAverage
   end with
next

t = 0
dt = 5e-15 's
while not multikey(1)
   'mouse handler
  mouseEvent = handleMouse(mouse)
  if (mouse.buttons <> -1) then
    if (mouseEvent = MOUSE_LB_PRESSED) then mouseDrag = 1
    if (mouseEvent = MOUSE_LB_RELEASED) then mouseDrag = 0
    if (mouseEvent = MOUSE_WHEEl_UP) then sg.scale *= 1.1
    if (mouseEvent = MOUSE_WHEEl_DOWN) then sg.scale /= 1.1
  end if

  if (mouseDrag) then
    sg.offset.x += mouse.posChange.x
    sg.offset.y -= mouse.posChange.y
  end if

   'drawing
   if (drawCounter = 10) then
      drawCounter = 0
      screenlock
      line (0,0)-(scrn_w-1, scrn_h-1),&h00000000,bf
      'sg.line_xy(0, -8e-9, 0, +8e-9, &h007f7f7f)
      'sg.line_xy(-8e-9, 0, +8e-9, 0, &h007f7f7f)
      for i = 0 to nSquares-1
         drawQuad(square(i), &h00ffff00)
      next
      for i = 0 to nBalls-1
         'sg.line(ball(i).p, vector_add(ball(i).p, vector_mul(ball(i).v, 0.005)), &h00ffffff)
         'sg.line(ball(i).p, vector_add(ball(i).p, vector_mul(ball(i).F, 2)), &h00ffffff)
         sg.circle(ball(i).p, ball(i).r, &h00ff00ff)
      next
      locate 1,1: print "Ekin  [J]"; Ekin;
      locate 2,1: print "Eelas [J]"; Eelas;
      locate 3,1: print "Etot  [J]"; Etot;
      screenunlock

'      if (colPairList.nItems > 0) then
'         print
'         print
'         colPairList.printList()
'         sleep
'      end if
      sleep 1,1
   end if
   drawCounter += 1

  'reset forces & Energy
   Ekin = 0
   Eelas = 0
   for i = 0 to nBalls-1
      ball(i).F.x = 0
      ball(i).F.y = 0
   next

   'calculate collision forces with squares
   for i = 0 to nBalls-1
      xBall = ball(i).p.x
      yBall = ball(i).p.y
      rBall = ball(i).r
      for k = 0 to nSquares-1
         collisionCount = 0
         for j = 0 to 3
            'A = square(k).line(j).A
            'B = square(k).line((j+1) and 3).A
            'result = closestPointOnLineSegment(ball(i).p, A, B, P)
            'result = cpols(ball(i).p, square(k), j, P)
            result = cpols(ball(i).p, square(k).line(j), P)
            'sg.circle(P, 188 * 1e-12, &h0000ffff)
            if (result = true) then 'on line segment
               dx = P.x - xBall
               dy = P.y - yBall
               dCenter = sqr(dx * dx + dy * dy)
               dEdge = dCenter - rBall
               if(dEdge < 0) then 'hitting a line
                  collisionCount += 1
                  F = kBall * (dEdge / dCenter)
                  ball(i).F.x += F * dx
                  ball(i).F.y += F * dy
                  Eelas += 0.5 * kBall * dEdge * dEdge '0.5*k*x^2
               end if
            end if
         next
         if (collisionCount = 0) then 'then maybe on corner?
            'loop all corners
            for j = 0 to 3
               dx = square(k).line(j).A.x - xBall
               dy = square(k).line(j).A.y - yBall
               dCenter = sqr(dx * dx + dy * dy)
               dEdge = dCenter - rBall
               if(dEdge < 0) then
                  F = kBall * (dEdge / dCenter)
                  ball(i).F.x += F * dx
                  ball(i).F.y += F * dy
                  Eelas += 0.5 * kBall * dEdge * dEdge '0.5*k*x^2
               end if
            next
         end if
      next
   next

   if (bruteForce = 0) then
   
      xSectDiv = square1size / nSectx
      ySectDiv = square1size / nSecty

      for yi = 0 to nSecty-1
         for xi = 0 to nSectx-1
            section(xi, yi).nItems = 0
         next
      next
      for i = 0 to nBalls-1
         with ball(i)

            xi1 = int(((.p.x - .r) + square1size / 2) / xSectDiv) 'left
            if (xi1 < 0) then xi1 = 0
            if (xi1 >= nSectx) then xi1 = nSectx-1

            xi2 = int(((.p.x + .r) + square1size / 2) / xSectDiv) 'right
            if (xi2 < 0) then xi2 = 0
            if (xi2 >= nSectx) then xi2 = nSectx-1

            yi1 = int(((.p.y - .r) + square1size / 2) / ySectDiv) 'left
            if (yi1 < 0) then yi1 = 0
            if (yi1 >= nSecty) then yi1 = nSecty-1

            yi2 = int(((.p.y + .r) + square1size / 2) / ySectDiv) 'right
            if (yi2 < 0) then yi2 = 0
            if (yi2 >= nSecty) then yi2 = nSecty-1

            section(xi1, yi1).addBall(i)
            if (xi1 = xi2) then
               if (yi1 <> yi2) then
                  section(xi1, yi2).addBall(i)
               end if
            else
               section(xi2, yi1).addBall(i)
               if (yi1 <> yi2) then
                  section(xi1, yi2).addBall(i)
                  section(xi2, yi2).addBall(i)
               end if
            end if

         end with
      next
            
      colPairList.nItems = 0
      'calculate collision forces between balls
      for yi = 0 to nSecty-1
         for xi = 0 to nSectx-1
            nList = section(xi, yi).nItems - 1 'number of balls in section
            
            for i2 = 0 to nList 'loop section list
               i = section(xi, yi).index(i2) 'get ball index
               xBall = ball(i).p.x
               yBall = ball(i).p.y
               rBall = ball(i).r
               for j2 = i2+1 to nList 'loop section list too
                  j = section(xi, yi).index(j2) 'get ball index too
                  dx = xBall - ball(j).p.x
                  dy = yBall - ball(j).p.y
                  dCenter = sqr(dx * dx + dy * dy)
                  dEdge = dCenter - (rBall + ball(j).r)
                  if(dEdge < 0) then 'BOOM!
                     'j > i Always!?! yes
                     if (colPairList.addPairIfNew(i, j) = true) then
                        F = kBall * (dEdge / dCenter)
                        ball(i).F.x -= F * dx
                        ball(i).F.y -= F * dy
                        ball(j).F.x += F * dx
                        ball(j).F.y += F * dy
                        Eelas += 0.5 * kBall * dEdge * dEdge '0.5*k*x^2
                     else
                        'pair handled already, should not collide twice
                     end if
                  end if
               next
            next
            
         next
      next

   else 'bruteForce = 1

      for i = 0 to nBalls-1
         xBall = ball(i).p.x
         yBall = ball(i).p.y
         rBall = ball(i).r
         for j = i+1 to nBalls-1
            dx = xBall - ball(j).p.x
            dy = yBall - ball(j).p.y
            dCenter = sqr(dx * dx + dy * dy)
            dEdge = dCenter - (rBall + ball(j).r)
            if(dEdge < 0) then
               F = kBall * (dEdge / dCenter)
               ball(i).F.x -= F * dx
               ball(i).F.y -= F * dy
               ball(j).F.x += F * dx
               ball(j).F.y += F * dy
               Eelas += 0.5 * kBall * dEdge * dEdge '0.5*k*x^2
            end if
         next
      next
      
   end if

   'calculate new positions
   for i = 0 to nBalls-1
      with ball(i)
         .a.x = (ball(i).F.x / ball(i).m)
         .a.y = (ball(i).F.y / ball(i).m) - g
         .v.x += .a.x * dt
         .v.y += .a.y * dt
         .p.x += .v.x * dt
         .p.y += .v.y * dt
         Ekin += 0.5 * ball(i).m * (.v.x * .v.x + .v.y * .v.y) 'E=m*v^2
      end with
   next
   Etot = Ekin + Eelas
      
   'ball(0).p.x = (mouse.pos.x - sg.offset.x) / sg.scale
   'ball(0).p.y = ((scrn_h - mouse.pos.y) - sg.offset.y) / sg.scale

   t += dt
   loopCounter += 1
   'if (loopCounter > 20000) then exit while

wend
print timer - startTime
print
for i = 0 to nTimers-1
   print tt(i).name, tt(i).total
next

sleep
end
deltarho[1859]
Posts: 2616
Joined: Jan 02, 2017 0:34
Location: UK

Re: Corona virus simulator

Postby deltarho[1859] » Mar 27, 2020 6:37

UEZ wrote:well done but the move of the pixels are too turbulent imho...

@ShawnLG

Have a look at 'Randomize , 5'. Not only do we get a 'slowing down' but we also get top drawer randomness.
UEZ
Posts: 627
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Postby UEZ » Mar 27, 2020 7:22

deltarho[1859] wrote:
UEZ wrote:well done but the move of the pixels are too turbulent imho...

@ShawnLG

Have a look at 'Randomize , 5'. Not only do we get a 'slowing down' but we also get top drawer randomness.

I would rather use a 2d Noise function to move the pixel if static moving vectors are not wanted.
UEZ
Posts: 627
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Postby UEZ » Mar 27, 2020 9:36

badidea wrote:
UEZ wrote:One of my first ideas was to create zones and check infection there only to speed it up but I don't know whether it was a good idea...

Yes, you can do a 'pre-selection'. You know that a particle or person in de top-right corner can never interact with bottom-left. Spilt the map in a grid. Only check within a grid-section, and neighbouring grid-sections. Make it flexible, so you can test what grid size is optimal. I did this before, I'll look up the code...

I found the code I wrote 8 years ago. A particle collision simulation. The variable bruteForce determines whether to check all particle against each other (1), or with pre-sectioning in a 8 x 8 grid (0). Will take me some time to understand my own code. I should have put more comments in. There is a nice feature: You can zoom in and out with the mouse-wheel and drag the view with right-mouse-button pressed.'


I increased the zones to 256 (2^8) and doubled the pixel size. Still over 30 fps when 90% are infected! See my 1st post in this topic for the code.
Be careful when changing the amount of zones! E.g. 128 zones will not cover the whole screen! Maybe I need to fix that part of the code...

It looks very nice. I will have a close look to your code..
badidea
Posts: 2149
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Mar 28, 2020 0:31

UEZ wrote:It looks very nice. I will have a close look to your code..

Don't look too close. I have my doubts about the efficiency of the grid-sectioning as I implemented it.
I assign each ball to a grid position, but when it overlaps a grid line or edge, i assign the ball to the neighbour sections as well. Then later on I keep a list of processed collisions (two prevent double collision). But when I think of it now, it is overly complicated. No need to assign ball to multiple section and keep a collision list. I could be wrong thou. I'll try a rewrite this weekend...
In the mean time, read this: https://gameprogrammingpatterns.com/spa ... m's-length
dodicat
Posts: 6692
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Corona virus simulator

Postby dodicat » Mar 28, 2020 10:17

I had implemented a line segment distance via cross products and all the vector operators +,-,*(a scaler), e.t.c. it worked OK also it worked for 3D.
But I replaced it with a much faster 2D version using only school co-ordinate geometry.
Also the fbgfx filled circle is not very good for small circles (line 226), so I replaced it with a custom circle.

Code: Select all

Screen 19,32
Dim Shared As Integer xres,yres
Dim Shared As Any Ptr row:row=Screenptr
Dim Shared As Integer pitch
Screeninfo xres,yres,,,pitch
Dim Shared As Long reds,greens
Type V2
    As Single x,y,dx,dy
    As Long  radius
    As Ulong c 'colour
    As Ulong m 'mass
    As Long recovery
End Type

Type Line2D
    As Single v1x,v1y,v2x,v2y
End Type

Function segment_distance(lx1 As Single, _
    ly1 As Single, _
    lx2 As Single, _
    ly2 As Single, _
    px As Single,_
    py As Single, _
    Byref ox As Single=0,_
    Byref oy As Single=0) As Single
    Dim As Single M1,M2,C1,C2,B
    B=(Lx2-Lx1):If B=0 Then B=1e-20
    M2=(Ly2-Ly1)/B:If M2=0 Then M2=1e-20
    M1=-1/M2
    C1=py-M1*px
    C2=(Ly1*Lx2-Lx1*Ly2)/B
    Var L1=((px-lx1)*(px-lx1)+(py-ly1)*(py-ly1)),L2=((px-lx2)*(px-lx2)+(py-ly2)*(py-ly2))
    Var a=((lx1-lx2)*(lx1-lx2) + (ly1-ly2)*(ly1-ly2))
    Var a1=a+L1
    Var a2=a+L2
    Var f1=a1>L2,f2=a2>L1
    If f1 Xor f2 Then
        Var d1=((px-Lx1)*(px-Lx1)+(py-Ly1)*(py-Ly1))
        Var d2=((px-Lx2)*(px-Lx2)+(py-Ly2)*(py-Ly2))
        If d1<d2 Then Ox=Lx1:Oy=Ly1 : Return Sqr(d1) Else  Ox=Lx2:Oy=Ly2:Return Sqr(d2)
    End If
    Var M=M1-M2:If M=0 Then M=1e-20
    Ox=(C2-C1)/(M1-M2)
    Oy=(M1*C2-M2*C1)/M
    Return Sqr((px-Ox)*(px-Ox)+(py-Oy)*(py-Oy))
End Function

'optimize detection to save cpu.
Function DetectBallCollisions(Byref b1 As V2,b2 As V2) As Single
    Dim As Single xdiff = b2.x-b1.x
    Dim As Single ydiff = b2.y-b1.y
    If Abs(xdiff) > b2.radius*2 Then Return 0
    If Abs(ydiff) > b2.radius*2 Then Return 0
    Var L=Sqr(xdiff*xdiff+ydiff*ydiff)
    If L<=(b2.radius+b1.radius) Then Function=L
End Function

Sub Check_BallCollisions(b() As v2)
    For n1 As Long=Lbound(b) To Ubound(b)-1
        For n2 As Long=n1+1 To Ubound(b)
            Dim As Single  L= DetectBallCollisions(b(n1),b(n2))
            If L Then
                Dim As Single  impulsex=(b(n1).x-b(n2).x)
                Dim As Single  impulsey=(b(n1).y-b(n2).y)
                Dim As Single ln=Sqr(impulsex*impulsex+impulsey*impulsey)
                impulsex/=ln'normalize the impulse
                impulsey/=ln
                'set one ball to nearest non overlap position
                b(n1).x=b(n2).x+(b(n2).radius+b(n1).radius)*impulsex
                b(n1).y=b(n2).y+(b(n2).radius+b(n1).radius)*impulsey
               
                Dim As Single  impactx=b(n1).dx-b(n2).dx
                Dim As Single  impacty=b(n1).dy-b(n2).dy
                Dim As Single  dot=impactx*impulsex+impacty*impulsey
                'handle mass
                Dim As Single  mn2=b(n1).m/(b(n1).m+b(n2).m),mn1=b(n2).m/(b(n1).m+b(n2).m)
                b(n1).dx-=dot*impulsex*2*mn1
                b(n1).dy-=dot*impulsey*2*mn1
                b(n2).dx+=dot*impulsex*2*mn2
                b(n2).dy+=dot*impulsey*2*mn2
                '=======  collisionds done =====
                '=====  colours only  ==========
                Static As Ulong white=Rgb(255,255,255)
                If b(n1).c=Rgb(200,0,0) Andalso Rnd>.5 Then b(n2).c=Rgb(200,0,0)
                b(n1).recovery+=1
                If b(n1).recovery>20 And n1<>1 And b(n1).c<>white Then  b(n1).recovery=0:b(n1).c=Rgb(0,200,0)
            End If
           
        Next n2
    Next n1
End Sub

Sub check_ball_to_line_collisions(LN() As Line2D, ball() As V2)
    For z As Integer=Lbound(ball) To Ubound(ball)
        For z2 As Integer=Lbound(Ln) To Ubound(Ln)
            Dim As V2 closepoint
            Var seperation=segment_distance(Ln(z2).v1x,Ln(z2).v1y,Ln(z2).v2x,Ln(z2).v2y,ball(z).x,ball(z).y,closepoint.x,closepoint.y)
            If seperation<=ball(z).radius Then
                Var impactx=-ball(z).dx
                Var impacty=-ball(z).dy
                Var impulsex=(closepoint.x-ball(z).x)/seperation
                Var impulsey=(closepoint.y-ball(z).y)/seperation
                ball(z).x=closepoint.x-ball(z).radius*impulsex
                ball(z).y=closepoint.y-ball(z).radius*impulsey
                Var dv=impactx*impulsex+impacty*impulsey
                ball(z).dx+= 2*dv*impulsex
                ball(z).dy+= 2*dv*impulsey
            End If
        Next z2
    Next z
End Sub

Sub drawline(L As Line2D,col As Ulong)
    Line(L.v1x,L.v1y)-(L.v2x,L.v2y),col
End Sub

Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
    Static As Double timervalue,_lastsleeptime,t3,frames
    Var t=Timer
    frames+=1
    If (t-t3)>=1 Then t3=t:fps=frames:frames=0
    Var sleeptime=_lastsleeptime+((1/myfps)-T+timervalue)*1000
    If sleeptime<1 Then sleeptime=1
    _lastsleeptime=sleeptime
    timervalue=T
    Return sleeptime
End Function

Function inpolygon(p1() As v2,Byval p2 As v2) As Integer
    #define Winder(L1,L2,p) ((L1.x-L2.x)*(p.y-L2.y)-(p.x-L2.x)*(L1.y-L2.y))
    Dim As Long index,nextindex,k=Ubound(p1)+1,wn
    For n As Long=1 To Ubound(p1)
        index=n Mod k:nextindex=(n+1) Mod k
        If nextindex=0 Then nextindex=1
        If p1(index).y<=p2.y Then
            If p1(nextindex).y>p2.y Andalso  Winder(p1(index),p1(nextindex),p2)>0 Then wn+=1
        Else
            If p1(nextindex).y<=p2.y Andalso Winder(p1(index),p1(nextindex),p2)<0 Then wn-=1
        End If
    Next n
    Return wn
End Function

Sub setuplines(segs() As Line2D,nodes() As v2)
    Dim As V2 Tmp(1 To 8)
    Redim  segs(1 To 7)
    For n As Integer=1 To 8
        Read Tmp(n).x
    Next n
    For n As Integer=1 To 8
        Read Tmp(n).y
    Next n
    For n As Integer=1 To 7
        segs(n)=Type<Line2D>(Tmp(n).x,Tmp(n).y,Tmp(n+1).x,Tmp(n+1).y)
    Next n
    For n As Long=1 To 7
        Redim Preserve nodes(1 To n)
        nodes(n)=Type(Tmp(n).x,tmp(n).y)
    Next
    'screen edge
    Redim Preserve segs(Lbound(segs) To Ubound(segs)+4)
    segs(8)=Type<Line2D>(30,30,xres-10,3)
    segs(9)=Type<Line2D>(xres-10,3,xres-10,yres-10)
    segs(10)=Type<Line2D>(xres-10,yres-10,10,yres-10)
    segs(11)=Type<Line2D>(10,yres-10,30,30)
End Sub

Sub setupballs(b() As v2,segs() As Line2D,nodes() As v2)
    Dim As Long flag,condition,ct
    For x As Long=40 To xres-40 Step 30
        For y As Long=40 To yres-40 Step 30
            ct+=1
            Redim Preserve b(1 To ct)
            b(1).c=Rgb(200,0,0)
            b(ct).x=x:b(ct).y=y
            b(ct).radius=4+Rnd*3
            b(ct).m=b(ct).radius^2
            Do
                b(ct).dx=(Rnd-Rnd)/2
                b(ct).dy=(Rnd-Rnd)/2
            Loop Until Abs(b(ct).dx)>.1 And Abs(b(ct).dy)>.1
            If inpolygon(nodes(),Type(x,y)) Then b(ct).c=Rgb(255,255,255) Else If b(ct).c<>Rgb(200,0,0) Then b(ct).c=Rgb(0,200,0)   
        Next y
    Next x
End Sub

Sub moveballs(b() As v2,Byref e As Long)
    e=0
    For z As Long=1 To Ubound(b)
        b(z).x+=b(z).dx
        b(z).y+=b(z).dy
        e+=.5*b(z).m*(b(z).dx*b(z).dx + b(z).dy*b(z).dy)
    Next z
End Sub

Sub drawlines(segs() As Line2D)
    For n As Long=Lbound(segs) To Ubound(segs)
        drawline(segs(n),Rgb(200,200,200))
    Next n
    Paint(1,1),Rgb(0,100,255),Rgb(200,200,200)
End Sub

Sub _circle(b As v2)
    #define incircle(cx,cy,radius,x,y) (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
    #define onscreen x>=0 And x<xres And y>.0 And y<yres
    #define putpixel(_x,_y,colour)    *Cptr(Ulong Ptr,row+ (_y)*pitch+ (_x) Shl 2)  =(colour)
    Dim As Ulong tc
    For x As Long=b.x-b.radius To b.x+b.radius
        For y As Long=b.y-b.radius To b.y+b.radius
            If incircle(b.x,b.y,b.radius,x,y) Andalso onscreen Then
                putpixel(x,y,b.c)
            End If
        Next
    Next
End Sub

Sub drawballs(b() As v2)
    reds=0
    greens=0
    Static As Ulong red=Rgb(200,0,0)
    Static As Ulong green=Rgb(0,200,0)
    For z As Integer=1 To Ubound(b)
        _circle(b(z))
        'Circle(b(z).x,b(z).y),b(z).radius,b(z).c,,,,f
        If b(z).c=red Then reds+=1
        If b(z).c=green Then greens+=1
    Next z
End Sub

Function start As Long
    Redim As Line2D segs()
    Redim As v2 nodes()
    setuplines(segs(),nodes())
    Redim As v2 b()
    setupballs(b(),segs(),nodes())
   
    Dim As Long fps,energy
   
    Do
        check_ball_to_line_collisions(segs(),b())
        Check_BallCollisions(b())
        moveballs(b(),energy)
        Screenlock
        Cls
        drawlines(segs())
        drawballs(b())
        Draw String(30,40),"FPS = " &fps
        Draw String(30,60),"Energy " &energy
       
        Draw String(30,80),"reds    " &reds
        Draw String(30,100),"greens " &greens
        Screenunlock
        Sleep regulate(100,fps),1
    Loop Until Len(Inkey)
    Sleep
    Return 0
End Function

End start

'the line segments
X_values:

Data _
247, 420, 654, 599, 470, 323, 210, 247

Y_values:

Data _
132, 78, 78, 417, 443, 451, 342, 132

 

The simulator bit is only a crude extra.
Tourist Trap
Posts: 2933
Joined: Jun 02, 2015 16:24

Re: Corona virus simulator

Postby Tourist Trap » Mar 28, 2020 11:33

dodicat wrote:I had implemented a line segment distance via cross products and all the vector operators +,-,*(a scaler), e.t.c. it worked OK also it worked for 3D.

Very nice.
I guess the immunized ones under heavy fences are the citizen of the crown of britain :)
UEZ
Posts: 627
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Postby UEZ » Mar 28, 2020 15:01

Very cool @dodicat.

@badidea: thanks for the link.
BasicScience
Posts: 489
Joined: Apr 18, 2008 4:09
Location: Los Angeles, CA
Contact:

Re: Corona virus simulator

Postby BasicScience » Mar 28, 2020 17:31

Nice work. Could be interesting to add two new features:

1) Motion is restricted to a maximum distance from initial position. This is a crude simulation of "stay home".

2) Simulate imperfect PPE. So each collision with an infected person has a defined probability of infecting the healthy one.
badidea
Posts: 2149
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Mar 29, 2020 0:06

While working on my own virus-outbreak-simulator, I got this hypnotizing graphics show.

Code: Select all

#Include "string.bi"

'===============================================================================

Type int2d
   As Integer x, y
   Declare Constructor
   Declare Constructor(x As Integer, y As Integer)
   Declare Operator Cast () As String
End Type

Constructor int2d
End Constructor

Constructor int2d(x As Integer, y As Integer)
   This.x = x : This.y = y
End Constructor

' "x, y"
Operator int2d.cast () As String
  Return Str(x) & "," & Str(y)
End Operator

Operator = (a As int2d, b As int2d) As boolean
   If a.x <> b.x Then Return false
   If a.y <> b.y Then Return false
   Return true
End Operator

Operator <> (a As int2d, b As int2d) As boolean
   If a.x = b.x And a.y = b.y Then Return false
   Return true
End Operator

' a + b
Operator + (a As int2d, b As int2d) As int2d
   Return Type(a.x + b.x, a.y + b.y)
End Operator

' a - b
Operator - (a As int2d, b As int2d) As int2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As int2d) As int2d
   Return Type(-a.x, -a.y)
End Operator

' a * b
Operator * (a As int2d, b As int2d) As int2d
   Return Type(a.x * b.x, a.y * b.y)
End Operator

' a * mul
Operator * (a As int2d, mul As Integer) As int2d
   Return Type(a.x * mul, a.y * mul)
End Operator

' a \ b
Operator \ (a As int2d, b As int2d) As int2d
   Return Type(a.x \ b.x, a.y \ b.y)
End Operator

' a \ div
Operator \ (a As int2d, divider As Integer) As int2d
   Return Type(a.x \ divider, a.y \ divider)
End Operator

'===============================================================================

Type sgl2d
   As Single x, y
   Declare Constructor
   Declare Constructor(x As Single, y As Single)
   Declare Operator Cast() As String
   Declare Function cross(b As sgl2d) As Single
   Declare Function lengthSqrd() As Single
   Declare Function dist(b As sgl2d) As Single
   Declare Function distSqrd(b As sgl2d) As Single
   Declare Function normalise() As sgl2d
End Type

Constructor sgl2d
End Constructor

Constructor sgl2d(x As Single, y As Single)
   This.x = x : This.y = y
End Constructor

Function sgl2d.cross(b As sgl2d) As Single
   Return This.x * b.y - This.y * b.x
End Function

Function sgl2d.lengthSqrd() As Single
   Return (This.x * This.x) + (This.y * This.y)
End Function

Function sgl2d.dist(b As sgl2d) As Single
   Dim As Single dx = This.x - b.x
   Dim As Single dy = This.y - b.y
   Return Sqr((dx * dx) + (dy * dy))
End Function

Function sgl2d.distSqrd(b As sgl2d) As Single
   Dim As Single dx = This.x - b.x
   Dim As Single dy = This.y - b.y
   Return (dx * dx) + (dy * dy)
End Function

Function sgl2d.normalise() As sgl2d
   Dim As Single length = Sqr((This.x * This.x) + (This.y * This.y))
   Return sgl2d(This.x / length, This.y / length)
End Function

' "x, y"
Operator sgl2d.cast() As String
   Return Str(x) & "," & Str(y)
End Operator

'---- operators ---

' distance / lenth
Operator Len (a As sgl2d) As Single
   Return Sqr(a.x * a.x + a.y * a.y)
End Operator

' a = b ?
Operator = (a As sgl2d, b As sgl2d) As boolean
   If a.x <> b.x Then Return false
   If a.y <> b.y Then Return false
   Return true
End Operator

' a != b ?
Operator <> (a As sgl2d, b As sgl2d) As boolean
   If a.x = b.x And a.y = b.y Then Return false
   Return true
End Operator

' a + b
Operator + (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x + b.x, a.y + b.y)
End Operator

' a - b
Operator - (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As sgl2d) As sgl2d
   Return Type(-a.x, -a.y)
End Operator

' a * b
Operator * (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x * b.x, a.y * b.y)
End Operator

' a * mul
Operator * (a As sgl2d, mul As Single) As sgl2d
   Return Type(a.x * mul, a.y * mul)
End Operator

' a / div
Operator / (a As sgl2d, div As Single) As sgl2d
   Return Type(a.x / div, a.y / div)
End Operator

'===============================================================================

'#include "vector.bi"
#Define MOUSE_IDLE 0
#Define MOUSE_POS_CHANGED 1
#Define MOUSE_LB_PRESSED 2
#Define MOUSE_LB_RELEASED 3
#Define MOUSE_RB_PRESSED 4
#Define MOUSE_RB_RELEASED 5
#Define MOUSE_MB_PRESSED 6
#Define MOUSE_MB_RELEASED 7
#Define MOUSE_WHEEL_UP 8
#Define MOUSE_WHEEL_DOWN 9

Type mouse_type
   Pos As int2d
   posChange As int2d
   wheel As Integer
   buttons As Integer
   lb As Integer 'left button
   rb As Integer 'right button
   mb As Integer 'middle button
End Type

Function handleMouse(ByRef mouse As mouse_type) As Integer
   Static previous As mouse_type
   Dim As Integer change = MOUSE_IDLE
   GetMouse mouse.pos.x, mouse.pos.y, mouse.wheel, mouse.buttons
   If (mouse.buttons = -1) Then
      mouse.lb = 0
      mouse.rb = 0
      mouse.mb = 0
      mouse.posChange.x = 0
      mouse.posChange.y = 0
   Else
      mouse.lb = (mouse.buttons And 1)
      mouse.rb = (mouse.buttons ShR 1) And 1
      mouse.mb = (mouse.buttons ShR 2) And 1
      If (previous.pos.x <> mouse.pos.x Or previous.pos.y <> mouse.pos.y) Then
         change = MOUSE_POS_CHANGED
      End If
      mouse.posChange.x = mouse.pos.x - previous.pos.x
      mouse.posChange.y = mouse.pos.y - previous.pos.y
      If (previous.buttons <> mouse.buttons) Then
         If (previous.lb = 0 And mouse.lb = 1) Then change = MOUSE_LB_PRESSED
         If (previous.lb = 1 And mouse.lb = 0) Then change = MOUSE_LB_RELEASED
         If (previous.rb = 0 And mouse.rb = 1) Then change = MOUSE_RB_PRESSED
         If (previous.rb = 1 And mouse.rb = 0) Then change = MOUSE_RB_RELEASED
         If (previous.mb = 0 And mouse.mb = 1) Then change = MOUSE_MB_PRESSED
         If (previous.mb = 1 And mouse.mb = 0) Then change = MOUSE_MB_RELEASED
      End If
      If (mouse.wheel > previous.wheel) Then change = MOUSE_WHEEl_UP
      If (mouse.wheel < previous.wheel) Then change = MOUSE_WHEEl_DOWN
      previous = mouse
   End If
   Return change
End Function

'===============================================================================

'#include "scaled_graphics.bi"
'Note: y+ = up, x+ = right, (0,0) = center
Type scaled_graphics_type
   Dim As Single scale = 1 ' = 1 / pixel_size 'pixels / meter
   'dim as int2d offset' = (scrn_w \ 2, h \ 2) 'offset in pixels
   Dim As sgl2d offset
   Dim As Integer w = -1, h = -1
   Dim As Integer wc = -1, hc = -1 'center x,y
   Declare Sub setScreen(w As Integer, h As Integer)
   Declare Sub setScaling(scale As Single, offset As sgl2d)
   Declare Sub clearScreen(c As ULong)
   Declare Function pos2screen(p As sgl2d) As int2d
   Declare Sub drawPixel(p As sgl2d, c As ULong)
   Declare Sub drawCircle(p As sgl2d, r As Single, c As ULong)
   Declare Sub drawCircleFilled(p As sgl2d, r As Single, c As ULong, cFill As ULong)
   Declare Sub drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)
   Declare Sub drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)
End Type

Sub scaled_graphics_type.setScreen(w As Integer, h As Integer)
   This.w = w 'width
   This.h = h 'height
   wc = w \ 2
   hc = h \ 2
   ScreenRes w, h, 32
   Width w \ 8, h \ 16 'bigger font
End Sub

Sub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)
   This.scale = scale
   This.offset = offset
End Sub

Sub scaled_graphics_type.clearScreen(c As ULong)
   Line(0, 0)-(w - 1, h - 1), c, bf
End Sub

Function scaled_graphics_type.pos2screen(p As sgl2d) As int2d
   Return int2d(Int(wc + (p.x - offset.x) * scale), h - Int(hc + (p.y - offset.y) * scale))
End Function

Sub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   PSet(posScrn.x, posScrn.y), c
End Sub

Sub scaled_graphics_type.drawCircle(p As sgl2d, r As Single, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   Circle(posScrn.x, posScrn.y), r * scale, c
End Sub

Sub scaled_graphics_type.drawCircleFilled(p As sgl2d, r As Single, c As ULong, cFill As ULong)
   Dim As int2d posScrn = pos2screen(p)
   Circle(posScrn.x, posScrn.y), r * scale, 0,,,,f
   Circle(posScrn.x, posScrn.y), r * scale, c
End Sub

Sub scaled_graphics_type.drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   Circle(posScrn.x, posScrn.y), r * scale, c, , , aspect
End Sub

Sub scaled_graphics_type.drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)
   Dim As int2d posScrn1 = pos2screen(p1)
   Dim As int2d posScrn2 = pos2screen(p2)
   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c
End Sub

'===============================================================================

'screen stuff
Const As Single PPM = 10.0 'pixels per meter
Const SCRN_W = 850, SCRN_H = 750

Dim Shared As scaled_graphics_type sg
sg.setScreen(SCRN_W, SCRN_H)
sg.setScaling(PPM, sgl2d(0, 0))

Type person_type
   p As sgl2d 'position [m]
   v As sgl2d 'velocity [m/s]
   r As Single 'radius [m]
End Type

Function drawUpdate(interval As Integer) As boolean
   Static As Integer counter = 0
   counter += 1
   If counter >= interval Then
      counter = 0
      Return TRUE
   Else
      Return FALSE
   End If
End Function

Function rndRangeSgl(min As Single, max As Single) As Single
   Return (Int((Rnd * (max - min) + min) * 5))/5
End Function

Const PERSONS = 4000
Const As Single MAP_X_MIN = -35, MAP_X_MAX = +35
Const As Single MAP_Y_MIN = -35, MAP_Y_MAX = +35

Dim As Integer mouseEvent, mouseDrag
Dim As mouse_type mouse
Dim As person_type person(0 To PERSONS-1)

'initialise persons
For i As Integer = 0 To UBound(person)
   person(i).p.x = 0 'rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm
   person(i).p.y = 0 'rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm
   person(i).v.x = rndRangeSgl(-5.0, +5.0) 'm/s
   person(i).v.y = rndRangeSgl(-5.0, +5.0) 'm/s
   person(i).r = 0.25
Next

Dim As Double tNow = Timer, tPrev = tNow, tStart = tNow
Dim As Double dt = tNow - tPrev
While Not MultiKey(1) 'escape
   mouseEvent = handleMouse(mouse)
      If (mouse.buttons <> -1) Then
      If (mouseEvent = MOUSE_LB_PRESSED) Then mouseDrag = 1
      If (mouseEvent = MOUSE_LB_RELEASED) Then mouseDrag = 0
      If (mouseEvent = MOUSE_WHEEl_UP) Then sg.scale *= 1.1
      If (mouseEvent = MOUSE_WHEEl_DOWN) Then sg.scale /= 1.1
   End If
   If (mouseDrag) Then
      sg.offset.x -= (mouse.posChange.x / PPM)
      sg.offset.y += (mouse.posChange.y / PPM)
   End If

   'update positions
   For i As Integer = 0 To UBound(person)
      person(i).p += person(i).v * dt
   Next
   'check wall collisions
   For i As Integer = 0 To UBound(person)
      If person(i).p.x < MAP_X_MIN Then person(i).v.x = +Abs(person(i).v.x)
      If person(i).p.x > MAP_X_MAX Then person(i).v.x = -abs(person(i).v.x)
      If person(i).p.y < MAP_Y_MIN Then person(i).v.y = +Abs(person(i).v.y)
      If person(i).p.y > MAP_Y_MAX Then person(i).v.y = -abs(person(i).v.y)
   Next
   
   If (drawUpdate(10) = TRUE) Then
      ScreenLock
      sg.clearScreen(0)
      For i As Integer = 0 To UBound(person)
         sg.drawCircle(person(i).p, person(i).r, &h00ff00ff)
      Next
      Locate 1,1: Print Format(tNow - tStart, "0.#0")
      ScreenUnLock
      Sleep 1,1 'don't sleep every loop
   End If
   tPrev = tNow
   tNow = Timer
   dt = tNow - tPrev
Wend
End

After exactly 350 seconds it return to initial state.

Return to “General”

Who is online

Users browsing this forum: No registered users and 4 guests