Corona virus simulator

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

Re: Corona virus simulator

Post by albert »

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: 408
Joined: Feb 01, 2007 16:54
Location: usa

Re: Corona virus simulator

Post by integer »

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

Re: Corona virus simulator

Post by UEZ »

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: Select all

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.[/quote]

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...

[quote="dodicat"]Not so much a simulation but a message.
[/quote]
Very nice idea. Keep well and fit.
ShawnLG
Posts: 142
Joined: Dec 25, 2008 20:21

Re: Corona virus simulator

Post by ShawnLG »

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: 862
Joined: May 05, 2015 5:35
Location: Germany

Re: Corona virus simulator

Post by grindstone »

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

Re: Corona virus simulator

Post by badidea »

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: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Corona virus simulator

Post by deltarho[1859] »

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: 972
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Post by UEZ »

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: 972
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Post by UEZ »

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: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Post by badidea »

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: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Corona virus simulator

Post by dodicat »

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: 2958
Joined: Jun 02, 2015 16:24

Re: Corona virus simulator

Post by Tourist Trap »

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: 972
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Post by UEZ »

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

Post by BasicScience »

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: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Post by badidea »

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.
Post Reply