Corona virus simulator

[Tourist Trap]

After exactly 350 seconds it return to initial state.

I can't say what exactly we are supposed to see. Impressive anyway :)

[UEZ]

@badidea: lol, looks like an infection of a borg collective.


@BasicScience: is there a special candidate you are asking?

[BasicScience]

@UEZ the changes I suggested could be added to any of the posted simulations. Would be interesting with yours, for example.

[BasicScience]

@Badidea - mesmerizing.

I see you are planning to implement a maximum radius of motion for each person. I guess the initial position should be a random location in the personal circle. Otherwise the motion vector will be along a radius line, and the reflection at the circle boarder will just be an oscillation back and forth along a diameter. I started looking into coding this but... calculating where a ray intersects the circle would require solving quadratic (computationally expensive). Of course this would be done only for those steps when the distance from origin > radius. Perhaps a reflection could be done at the end of this step just outside the circle. The boundary would be fuzzy (but real life is fuzzy too), and it would save having to solve a quadratic.

[BasicScience]

Pseudo-random movement within a circle for "safer at home" simulation

Code: Select all

Type __person
    dim as single x_init, y_init, x , y, sqradius
    dim as single  Vx, Vy
    dim as ulong   c
END TYPE

DIM as integer N_people = 2, iPsize = 5

ScreenRes 640, 480, 32

DIM Shared as __person person(N_people-1)
Randomize

'Initialize
Person(0).x_init = 100
person(0).y_init = 100
person(0).Vx = 2*(rnd-.5)
person(0).Vy = 2*(rnd-.5)
person(0).c      = &hFFFFFFFF
person(0).sqradius = 2500

person(1).x_init = 300
person(1).y_init = 200
Person(1).Vx = (rnd-.5)
Person(1).Vy = (rnd-.5)
Person(1).c      = &hFFFFFFFF
Person(1).sqradius = 6400

For i as integer = 0 to N_People -1
    Person(i).x = Person(i).x_init + sqr(person(i).sqradius)*(rnd-.5)
    person(i).y = person(i).y_init + sqr(person(i).sqradius)*(rnd-.5)
    Line (Person(i).x, Person(i).y)-(Person(i).x + iPSize, Person(i).y + iPSize), Person(i).c, BF 
    circle (person(i).x_init, person(i).y_init),sqr(person(i).sqradius), rgb(255,0,0)
Next i

do until inkey = "" : loop

DIM as integer cnt = 0

Do
    CLS
    
    For i as integer = 0 to N_people-1
        IF CNT > 200 then
            Person(i).Vx = Person(i).Vx * sgn(rnd-.5)
            Person(i).Vy = Person(i).Vy * sgn(rnd-.5)
            Cnt = 0
        End IF
        Cnt += 1
        
        Person(i).x += Person(i).Vx 
        Person(i).y += Person(i).Vy
        
        IF (Person(i).x - Person(i).x_init) * (person(i).x - person(i).x_init) + (person(i).y - Person(i).y_init) * (person(i).y -person(i).y_init) > person(i).sqradius then
            
            'Get back inside circle
            Person(i).x -= Person(i).Vx
            Person(i).y -= Person(i).Vy
            
            'Assume current x,y is on circle (not really true)
            'Computer new velocities base on reflection to tangent
            
            DIM as single Nx, Ny, scale  ' normal
            Nx = -(Person(i).x - Person(i).X_init)
            Ny = -(Person(i).y - Person(i).Y_init)

            scale = 2*(Person(i).Vx * Nx + Person(i).Vy * Ny)/(Nx*Nx + Ny*Ny)

            Person(i).Vx = Person(i).Vx - scale * Nx
            Person(i).Vy = Person(i).Vy - Scale * Ny
            
            Person(i).x += Person(i).Vx
            Person(i).y += Person(i).Vy

        End if
        Line (Person(i).x, Person(i).y)-(Person(i).x + iPSize, Person(i).y + iPSize), Person(i).c, BF 
        circle (person(i).x_init, person(i).y_init),sqr(person(i).sqradius), rgb(255,0,0)
    next i
    Sleep 5
Loop while inkey = ""

sleep

[badidea]

Virus outbreak form the patient zero perspective. <space> to trigger. Uses a semi-random walk.
To do:
* Sectioning to handle more persons
* Social distancing control

Code: Select all

#Include "fbgfx.bi"
#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

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

#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

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

'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 drawCircleFilled(p As sgl2d, r As Single, c 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)
	Dim As int2d posScrn = pos2screen(p)
	Circle(posScrn.x, posScrn.y), r * scale, c,,,,f
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

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

Const As Single INFEC_DIST = 0.20 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY 

Const PERSONS = 500
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX

Const DRAW_INTERVAL = 1 '1 to 10 
Const FOLLOW_PATIENT_0 = TRUE 'FALSE

'screen stuff
Const As Single PPM = 20.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.setScreen(SCRN_W, SCRN_H)
sg.setScaling(PPM, sgl2d(0, 0))

Enum HEALTH_STATE
	ST_INIT '0
	ST_INFEC '1
	ST_RECOV '2
	ST_DEAD '3
	ST_LAST 'number of states
End Enum

Dim As Integer statCounters(ST_LAST - 1)
Dim As ULong stateColor(ST_LAST - 1) = _
	{RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}

Type person_type
	p As sgl2d 'position [m]
	v As sgl2d 'velocity [m/s]
	r As Single 'radius [m]
	sickEndTime As Single 'not double!
	state As Integer 'health
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 Rnd() * (max - min) + min
End Function

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

Randomize Timer
'initialise persons
For i As Integer = 0 To UBound(person)
	person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm
	person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm
	person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).r = 0.25
	'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))
	person(i).state = ST_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov

While Not MultiKey(FB.SC_ESCAPE)
	'trigger outbreak
	If MultiKey(FB.SC_SPACE) Then
		person(0).state = ST_INFEC
		person(0).sickEndTime = simulTime _
			+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
	End If
	'zoom / drag view by mouse
	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
	'patient 0 perpective
	If FOLLOW_PATIENT_0 = TRUE Then
		sg.offset.x = (person(0).p.x)
		sg.offset.y = (person(0).p.y)
	End If
	'update positions
	For i As Integer = 0 To UBound(person)
		person(i).p += person(i).v * timeStep
	Next
	'random direction change
	For i As Integer = 1 To 5
		Dim As Integer iPerson = Int(rndRangeSgl(0, PERSONS))
		If Rnd() > 0.5 Then
			person(iPerson).v.x = -person(iPerson).v.x
		Else
			person(iPerson).v.y = -person(iPerson).v.y
		End If
	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
	'check end of sickness
	For i As Integer = 0 To UBound(person)
		If person(i).state = ST_INFEC Then
			If simulTime > person(i).sickEndTime Then
				If Rnd() < MORTALITY Then
					person(i).state = ST_DEAD
					person(i).v = sgl2d(0, 0)
				Else
					person(i).state = ST_RECOV
					'person(i).v /= SICK_SPEED_FACTOR
				End If
			End If
		End If
	Next
	'check for disease transmission (iSrc = source, iTar = target)
	For iSrc As Integer = 0 To UBound(person)
		If person(iSrc).state = ST_INFEC Then
			For iTar As Integer = 0 To UBound(person)
				If person(iTar).state = ST_INIT Then
					If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
						person(iTar).state = ST_INFEC
						person(iTar).sickEndTime = simulTime _
							+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
						'person(iTar).v *= SICK_SPEED_FACTOR
					End If
				End If
			Next
		End If
	Next
	'clear stats
	For i As Integer = 0 To UBound(statCounters)
		statCounters(i) = 0
	Next
	'update stats
	For i As Integer = 0 To UBound(person)
		statCounters(person(i).state) += 1
	Next
	'draw world
	If (drawUpdate(DRAW_INTERVAL) = TRUE) Then
		ScreenLock
		sg.clearScreen(RGB(250, 250, 250))
		sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))
		For i As Integer = 0 To UBound(person)
			Dim As ULong c = stateColor(person(i).state)
			sg.drawCircleFilled(person(i).p, person(i).r, c)
		Next
		Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)
		Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)
		Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)
		Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)
		Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)
		ScreenUnLock
		Sleep 1,1 'don't sleep every loop
	End If
	simulTime += timeStep
	simulDays = simulTime * DAY_PER_SEC
Wend
End

[albert]

I came up with a cure for Covid-19

You take the white blood cells from persons who have recovered the fastest..
And transplant them into persons who aren't recovering well.. ( people in the hospitals )

It would also work for the regular influenza that kills far more people each year..

[badidea]

A version with grid based spacial partitioning. I hope that I did not make an error, with 6 nested for-loops and a bunch of if-statements, a bit nasty.
It is a bit heavy on the drawing, so I only display 1 in 10 frames. Can be speed up further with macros, ditching the classes and testing the optimal spatial sectioning grid size.

Code: Select all

#Include "fbgfx.bi"
#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

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

'list can grow, but never shrink, for performance, non-sorted
Type grow_list
	Dim As Integer index(Any)
	Private:
	Dim As Integer numItems
	Public:
	Declare Constructor(startSize As Integer)
	Declare Constructor()
	Declare Destructor()
	Declare Sub Add(newIndex As Integer)
	Declare Sub empty()
	Declare Function numAlloc() As Integer
	Declare Function numInUse() As Integer
	Declare Sub show()
	'non-list methods
End Type

Constructor grow_list(startSize As Integer)
	If startSize > 0 Then
		ReDim index(startSize - 1)
	End If
End Constructor

Constructor grow_list()
	This.constructor(0)
End Constructor

Destructor grow_list()
	Erase(index)
End Destructor

Sub grow_list.add(newIndex As Integer)
	Dim As Integer ub = UBound(index)
	'if list is full, increase list size by 1
	If numItems = ub + 1 Then
		ReDim Preserve index(numItems)
	End If
	index(numItems) = newIndex
	numItems += 1
End Sub

Sub grow_list.empty()
	numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
	Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
	Return numItems
End Function

'for debugging
Sub grow_list.show()
	Print "--- " & numInUse() & " / " & numAlloc() & " ---"
	For i As Integer = 0 To numItems - 1
		Print i, index(i)
	Next
End Sub

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

#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

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

'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 drawCircleFilled(p As sgl2d, r As Single, c 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)
	Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)
	Declare Sub drawRectFilled(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)
	Dim As int2d posScrn = pos2screen(p)
	Circle(posScrn.x, posScrn.y), r * scale, c,,,,f
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

Sub scaled_graphics_type.drawRect(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, b
End Sub

Sub scaled_graphics_type.drawRectFilled(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, bf
End Sub

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

Const As Single INFEC_DIST = 0.25 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY 

Const PERSONS = 10000
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10 
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)
'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]

'screen stuff
Const As Single PPM = 10.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.setScreen(SCRN_W, SCRN_H)
sg.setScaling(PPM, sgl2d(0, 0))

Enum HEALTH_STATE
	ST_INIT '0
	ST_INFEC '1
	ST_RECOV '2
	ST_DEAD '3
	ST_LAST 'number of states
End Enum

Dim As Integer statCounters(ST_LAST - 1)
Dim As ULong stateColor(ST_LAST - 1) = _
	{RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}

Type person_type
	p As sgl2d 'position [m]
	v As sgl2d 'velocity [m/s]
	r As Single 'radius [m]
	sickEndTime As Single 'not double!
	state As Integer 'health
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 Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
	If value < min Then value = min
	If value > max Then value = max
End Sub

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

Randomize Timer
'initialise persons
For i As Integer = 0 To UBound(person)
	person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm
	person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm
	person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).r = 0.25
	'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))
	person(i).state = ST_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

While Not MultiKey(FB.SC_ESCAPE)
	'trigger outbreak
	If MultiKey(FB.SC_SPACE) Then
		person(0).state = ST_INFEC
		person(0).sickEndTime = simulTime _
			+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
	End If
	'zoom / drag view by mouse
	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
	'patient 0 perpective
	If FOLLOW_PATIENT_0 = TRUE Then
		sg.offset.x = (person(0).p.x)
		sg.offset.y = (person(0).p.y)
	End If
	'update positions
	For i As Integer = 0 To UBound(person)
		person(i).p += person(i).v * timeStep
	Next
	'random direction change
	'~ for i as integer = 1 to 5
		'~ dim as integer iPerson = int(rndRangeSgl(0, PERSONS))
		'~ if rnd() > 0.5 then
			'~ person(iPerson).v.x = -person(iPerson).v.x
		'~ else
			'~ person(iPerson).v.y = -person(iPerson).v.y
		'~ end if
	'~ 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
	'check end of sickness
	For i As Integer = 0 To UBound(person)
		If person(i).state = ST_INFEC Then
			If simulTime > person(i).sickEndTime Then
				If Rnd() < MORTALITY Then
					person(i).state = ST_DEAD
					person(i).v = sgl2d(0, 0)
				Else
					person(i).state = ST_RECOV
					'person(i).v /= SICK_SPEED_FACTOR
				End If
			End If
		End If
	Next
	'clear sectors
	For xi As Integer = 0 To NUM_SECT_X - 1
		For yi As Integer = 0 To NUM_SECT_Y - 1
			sector(xi, yi).empty()
		Next
	Next
	'assign persons to sectors
	For i As Integer = 0 To UBound(person)
		Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))
		Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))
		limitInt(xi, 0, NUM_SECT_X - 1)
		limitInt(yi, 0, NUM_SECT_Y - 1)
		sector(xi, yi).add(i)
	Next
	If sectioning = TRUE Then
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				'loop source in 1 sector
				For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1
					Dim As Integer iSrc = sector(xi, yi).index(iiSrc)
					If person(iSrc).state = ST_INFEC Then
						'check against targets in near sectors, including this sector
						For xdi As Integer = -1 To +1
							If xi + xdi < 0 Then Continue For
							If xi + xdi >= NUM_SECT_X Then Continue For
							For ydi As Integer = -1 To +1
								If yi + ydi < 0 Then Continue For
								If yi + ydi >= NUM_SECT_Y Then Continue For
								'loop targets with 1 (near) sector
								For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1
									Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)
									If person(iTar).state = ST_INIT Then
										If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
											person(iTar).state = ST_INFEC
											person(iTar).sickEndTime = simulTime _
												+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
											'person(iTar).v *= SICK_SPEED_FACTOR
										End If
									End If
								Next
							Next
						Next
					End If
				Next
			Next
		Next
	Else
		'case: sectioning = FALSE
		'check for disease transmission (iSrc = source, iTar = target)
		For iSrc As Integer = 0 To UBound(person)
			If person(iSrc).state = ST_INFEC Then
				For iTar As Integer = 0 To UBound(person)
					If person(iTar).state = ST_INIT Then
						If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
							person(iTar).state = ST_INFEC
							person(iTar).sickEndTime = simulTime _
								+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
							'person(iTar).v *= SICK_SPEED_FACTOR
						End If
					End If
				Next
			End If
		Next
	End If
	'clear stats
	For i As Integer = 0 To UBound(statCounters)
		statCounters(i) = 0
	Next
	'update stats
	For i As Integer = 0 To UBound(person)
		statCounters(person(i).state) += 1
	Next
	'draw world
	If (drawUpdate(DRAW_INTERVAL) = TRUE) Then
		ScreenLock
		sg.clearScreen(RGB(250, 250, 250))
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN
				Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN
				sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), RGB(0, 0, 0))
				Dim As int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))
				Draw String (scrnPos.x + 2, scrnPos.y - 16), Str(sector(xi, yi).numInUse), 0
				'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0
			Next
		Next
		sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))
		For i As Integer = 0 To UBound(person)
			Dim As ULong c = stateColor(person(i).state)
			sg.drawCircleFilled(person(i).p, person(i).r, c)
		Next
		Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)
		Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)
		Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)
		Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)
		Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)
		ScreenUnLock
		Sleep 1,1 'don't sleep every loop
	End If
	simulTime += timeStep
	simulDays = simulTime * DAY_PER_SEC
Wend
End

[UEZ]

A version with grid based spacial partitioning. I hope that I did not make an error, with 6 nested for-loops and a bunch of if-statements, a bit nasty.
It is a bit heavy on the drawing, so I only display 1 in 10 frames. Can be speed up further with macros, ditching the classes and testing the optimal spatial sectioning grid size.

Code: Select all

#Include "fbgfx.bi"
#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

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

'list can grow, but never shrink, for performance, non-sorted
Type grow_list
	Dim As Integer index(Any)
	Private:
	Dim As Integer numItems
	Public:
	Declare Constructor(startSize As Integer)
	Declare Constructor()
	Declare Destructor()
	Declare Sub Add(newIndex As Integer)
	Declare Sub empty()
	Declare Function numAlloc() As Integer
	Declare Function numInUse() As Integer
	Declare Sub show()
	'non-list methods
End Type

Constructor grow_list(startSize As Integer)
	If startSize > 0 Then
		ReDim index(startSize - 1)
	End If
End Constructor

Constructor grow_list()
	This.constructor(0)
End Constructor

Destructor grow_list()
	Erase(index)
End Destructor

Sub grow_list.add(newIndex As Integer)
	Dim As Integer ub = UBound(index)
	'if list is full, increase list size by 1
	If numItems = ub + 1 Then
		ReDim Preserve index(numItems)
	End If
	index(numItems) = newIndex
	numItems += 1
End Sub

Sub grow_list.empty()
	numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
	Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
	Return numItems
End Function

'for debugging
Sub grow_list.show()
	Print "--- " & numInUse() & " / " & numAlloc() & " ---"
	For i As Integer = 0 To numItems - 1
		Print i, index(i)
	Next
End Sub

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

#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

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

'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 drawCircleFilled(p As sgl2d, r As Single, c 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)
	Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)
	Declare Sub drawRectFilled(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)
	Dim As int2d posScrn = pos2screen(p)
	Circle(posScrn.x, posScrn.y), r * scale, c,,,,f
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

Sub scaled_graphics_type.drawRect(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, b
End Sub

Sub scaled_graphics_type.drawRectFilled(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, bf
End Sub

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

Const As Single INFEC_DIST = 0.25 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY 

Const PERSONS = 10000
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10 
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)
'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]

'screen stuff
Const As Single PPM = 10.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.setScreen(SCRN_W, SCRN_H)
sg.setScaling(PPM, sgl2d(0, 0))

Enum HEALTH_STATE
	ST_INIT '0
	ST_INFEC '1
	ST_RECOV '2
	ST_DEAD '3
	ST_LAST 'number of states
End Enum

Dim As Integer statCounters(ST_LAST - 1)
Dim As ULong stateColor(ST_LAST - 1) = _
	{RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}

Type person_type
	p As sgl2d 'position [m]
	v As sgl2d 'velocity [m/s]
	r As Single 'radius [m]
	sickEndTime As Single 'not double!
	state As Integer 'health
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 Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
	If value < min Then value = min
	If value > max Then value = max
End Sub

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

Randomize Timer
'initialise persons
For i As Integer = 0 To UBound(person)
	person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm
	person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm
	person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).r = 0.25
	'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))
	person(i).state = ST_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

While Not MultiKey(FB.SC_ESCAPE)
	'trigger outbreak
	If MultiKey(FB.SC_SPACE) Then
		person(0).state = ST_INFEC
		person(0).sickEndTime = simulTime _
			+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
	End If
	'zoom / drag view by mouse
	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
	'patient 0 perpective
	If FOLLOW_PATIENT_0 = TRUE Then
		sg.offset.x = (person(0).p.x)
		sg.offset.y = (person(0).p.y)
	End If
	'update positions
	For i As Integer = 0 To UBound(person)
		person(i).p += person(i).v * timeStep
	Next
	'random direction change
	'~ for i as integer = 1 to 5
		'~ dim as integer iPerson = int(rndRangeSgl(0, PERSONS))
		'~ if rnd() > 0.5 then
			'~ person(iPerson).v.x = -person(iPerson).v.x
		'~ else
			'~ person(iPerson).v.y = -person(iPerson).v.y
		'~ end if
	'~ 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
	'check end of sickness
	For i As Integer = 0 To UBound(person)
		If person(i).state = ST_INFEC Then
			If simulTime > person(i).sickEndTime Then
				If Rnd() < MORTALITY Then
					person(i).state = ST_DEAD
					person(i).v = sgl2d(0, 0)
				Else
					person(i).state = ST_RECOV
					'person(i).v /= SICK_SPEED_FACTOR
				End If
			End If
		End If
	Next
	'clear sectors
	For xi As Integer = 0 To NUM_SECT_X - 1
		For yi As Integer = 0 To NUM_SECT_Y - 1
			sector(xi, yi).empty()
		Next
	Next
	'assign persons to sectors
	For i As Integer = 0 To UBound(person)
		Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))
		Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))
		limitInt(xi, 0, NUM_SECT_X - 1)
		limitInt(yi, 0, NUM_SECT_Y - 1)
		sector(xi, yi).add(i)
	Next
	If sectioning = TRUE Then
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				'loop source in 1 sector
				For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1
					Dim As Integer iSrc = sector(xi, yi).index(iiSrc)
					If person(iSrc).state = ST_INFEC Then
						'check against targets in near sectors, including this sector
						For xdi As Integer = -1 To +1
							If xi + xdi < 0 Then Continue For
							If xi + xdi >= NUM_SECT_X Then Continue For
							For ydi As Integer = -1 To +1
								If yi + ydi < 0 Then Continue For
								If yi + ydi >= NUM_SECT_Y Then Continue For
								'loop targets with 1 (near) sector
								For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1
									Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)
									If person(iTar).state = ST_INIT Then
										If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
											person(iTar).state = ST_INFEC
											person(iTar).sickEndTime = simulTime _
												+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
											'person(iTar).v *= SICK_SPEED_FACTOR
										End If
									End If
								Next
							Next
						Next
					End If
				Next
			Next
		Next
	Else
		'case: sectioning = FALSE
		'check for disease transmission (iSrc = source, iTar = target)
		For iSrc As Integer = 0 To UBound(person)
			If person(iSrc).state = ST_INFEC Then
				For iTar As Integer = 0 To UBound(person)
					If person(iTar).state = ST_INIT Then
						If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
							person(iTar).state = ST_INFEC
							person(iTar).sickEndTime = simulTime _
								+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
							'person(iTar).v *= SICK_SPEED_FACTOR
						End If
					End If
				Next
			End If
		Next
	End If
	'clear stats
	For i As Integer = 0 To UBound(statCounters)
		statCounters(i) = 0
	Next
	'update stats
	For i As Integer = 0 To UBound(person)
		statCounters(person(i).state) += 1
	Next
	'draw world
	If (drawUpdate(DRAW_INTERVAL) = TRUE) Then
		ScreenLock
		sg.clearScreen(RGB(250, 250, 250))
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN
				Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN
				sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), RGB(0, 0, 0))
				Dim As int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))
				Draw String (scrnPos.x + 2, scrnPos.y - 16), Str(sector(xi, yi).numInUse), 0
				'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0
			Next
		Next
		sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))
		For i As Integer = 0 To UBound(person)
			Dim As ULong c = stateColor(person(i).state)
			sg.drawCircleFilled(person(i).p, person(i).r, c)
		Next
		Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)
		Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)
		Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)
		Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)
		Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)
		ScreenUnLock
		Sleep 1,1 'don't sleep every loop
	End If
	simulTime += timeStep
	simulDays = simulTime * DAY_PER_SEC
Wend
End

Pretty cool.

You and dodicat using often the operator within you code. That's apparently a good way to make things easier. I need to understand it first...

[badidea]

I am running a comparison With and Without spacial partitioning (sectors in code).

I have added accumulative timers for simulation steps (excluding rendering):
step 1 = irrelevant (mouse and keyboard handling)
step 2 = update positions
step 3 = wall collisions
step 4 = end of sickness check
step 5 = irrelevant (clear sectors)
step 6 = assign persons to sectors
step 7 = virus transmission check
step 8 = irrelevant (update statistics)
step 9 = not a step, sum of timers above

Result with with spacial partitioning:


Result with without spacial partitioning:
(I let you know in an hour or more, going very slow, now at day 4 after 1400 seconds)
Done, that took some time (step 5 and 6 were not needed, but not significant):


The time from first outbreak to last patient recovered (or dead) is the same. This make me confident the the 'fast code' is correct.

Code: Select all

#Include "fbgfx.bi"
#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

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

'list can grow, but never shrink, for performance, non-sorted
Type grow_list
	Dim As Integer index(Any)
	Private:
	Dim As Integer numItems
	Public:
	Declare Constructor(startSize As Integer)
	Declare Constructor()
	Declare Destructor()
	Declare Sub Add(newIndex As Integer)
	Declare Sub empty()
	Declare Function numAlloc() As Integer
	Declare Function numInUse() As Integer
	Declare Sub show()
	'non-list methods
End Type

Constructor grow_list(startSize As Integer)
	If startSize > 0 Then
		ReDim index(startSize - 1)
	End If
End Constructor

Constructor grow_list()
	This.constructor(0)
End Constructor

Destructor grow_list()
	Erase(index)
End Destructor

Sub grow_list.add(newIndex As Integer)
	Dim As Integer ub = UBound(index)
	'if list is full, increase list size by 1
	If numItems = ub + 1 Then
		ReDim Preserve index(numItems)
	End If
	index(numItems) = newIndex
	numItems += 1
End Sub

Sub grow_list.empty()
	numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
	Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
	Return numItems
End Function

'for debugging
Sub grow_list.show()
	Print "--- " & numInUse() & " / " & numAlloc() & " ---"
	For i As Integer = 0 To numItems - 1
		Print i, index(i)
	Next
End Sub

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

#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

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

'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 drawCircleFilled(p As sgl2d, r As Single, c 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)
	Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)
	Declare Sub drawRectFilled(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)
	Dim As int2d posScrn = pos2screen(p)
	Circle(posScrn.x, posScrn.y), r * scale, c,,,,f
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

Sub scaled_graphics_type.drawRect(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, b
End Sub

Sub scaled_graphics_type.drawRectFilled(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, bf
End Sub

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

Const As Single INFEC_DIST = 0.25 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY 

Const PERSONS = 10000
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10 
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)
'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]

'screen stuff
Const As Single PPM = 10.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.setScreen(SCRN_W, SCRN_H)
sg.setScaling(PPM, sgl2d(0, 0))

Enum HEALTH_STATE
	ST_INIT '0
	ST_INFEC '1
	ST_RECOV '2
	ST_DEAD '3
	ST_LAST 'number of states
End Enum

Dim As Integer statCounters(ST_LAST - 1)
Dim As ULong stateColor(ST_LAST - 1) = _
	{RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}

Type person_type
	p As sgl2d 'position [m]
	v As sgl2d 'velocity [m/s]
	r As Single 'radius [m]
	sickEndTime As Single 'not double!
	state As Integer 'health
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 Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
	If value < min Then value = min
	If value > max Then value = max
End Sub

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

Randomize 1234
'initialise persons
For i As Integer = 0 To UBound(person)
	person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm
	person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm
	person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).r = 0.25
	'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))
	person(i).state = ST_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

sg.offset.x = -15 'move view for number
Dim As Double tNow(0 To 9), tAcc(0 To 9)
Dim As Integer trigger = 1

While Not MultiKey(FB.SC_ESCAPE)
	tNow(0) = Timer

	'trigger outbreak
	If MultiKey(FB.SC_SPACE) Then trigger = 1
	If trigger = 1 Then
		trigger = 0
		person(0).state = ST_INFEC
		person(0).sickEndTime = simulTime _
			+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
	End If
	'zoom / drag view by mouse
	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
	'patient 0 perpective
	If FOLLOW_PATIENT_0 = TRUE Then
		sg.offset.x = (person(0).p.x)
		sg.offset.y = (person(0).p.y)
	End If

	tNow(1) = Timer

	'update positions
	For i As Integer = 0 To UBound(person)
		person(i).p += person(i).v * timeStep
	Next
	'random direction change
	'~ for i as integer = 1 to 5
		'~ dim as integer iPerson = int(rndRangeSgl(0, PERSONS))
		'~ if rnd() > 0.5 then
			'~ person(iPerson).v.x = -person(iPerson).v.x
		'~ else
			'~ person(iPerson).v.y = -person(iPerson).v.y
		'~ end if
	'~ next

	tNow(2) = Timer

	'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

	tNow(3) = Timer

	'check end of sickness
	For i As Integer = 0 To UBound(person)
		If person(i).state = ST_INFEC Then
			If simulTime > person(i).sickEndTime Then
				If Rnd() < MORTALITY Then
					person(i).state = ST_DEAD
					person(i).v = sgl2d(0, 0)
				Else
					person(i).state = ST_RECOV
					'person(i).v /= SICK_SPEED_FACTOR
				End If
			End If
		End If
	Next
	
	tNow(4) = Timer

	'clear sectors
	For xi As Integer = 0 To NUM_SECT_X - 1
		For yi As Integer = 0 To NUM_SECT_Y - 1
			sector(xi, yi).empty()
		Next
	Next
	
	tNow(5) = Timer

	'assign persons to sectors
	For i As Integer = 0 To UBound(person)
		Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))
		Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))
		limitInt(xi, 0, NUM_SECT_X - 1)
		limitInt(yi, 0, NUM_SECT_Y - 1)
		sector(xi, yi).add(i)
	Next

	tNow(6) = Timer

	'spread the virus
	If sectioning = TRUE Then
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				'loop source in 1 sector
				For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1
					Dim As Integer iSrc = sector(xi, yi).index(iiSrc)
					If person(iSrc).state = ST_INFEC Then
						'check against targets in near sectors, including this sector
						For xdi As Integer = -1 To +1
							If xi + xdi < 0 Then Continue For
							If xi + xdi >= NUM_SECT_X Then Continue For
							For ydi As Integer = -1 To +1
								If yi + ydi < 0 Then Continue For
								If yi + ydi >= NUM_SECT_Y Then Continue For
								'loop targets with 1 (near) sector
								For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1
									Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)
									If person(iTar).state = ST_INIT Then
										If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
											person(iTar).state = ST_INFEC
											person(iTar).sickEndTime = simulTime _
												+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
											'person(iTar).v *= SICK_SPEED_FACTOR
										End If
									End If
								Next
							Next
						Next
					End If
				Next
			Next
		Next
	Else
		'case: sectioning = FALSE
		'check for disease transmission (iSrc = source, iTar = target)
		For iSrc As Integer = 0 To UBound(person)
			If person(iSrc).state = ST_INFEC Then
				For iTar As Integer = 0 To UBound(person)
					If person(iTar).state = ST_INIT Then
						If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
							person(iTar).state = ST_INFEC
							person(iTar).sickEndTime = simulTime _
								+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
							'person(iTar).v *= SICK_SPEED_FACTOR
						End If
					End If
				Next
			End If
		Next
	End If

	tNow(7) = Timer

	'clear stats
	For i As Integer = 0 To UBound(statCounters)
		statCounters(i) = 0
	Next
	'update stats
	For i As Integer = 0 To UBound(person)
		statCounters(person(i).state) += 1
	Next
	
	tNow(8) = Timer

	'draw world
	If (drawUpdate(DRAW_INTERVAL) = TRUE) Then
		ScreenLock
		sg.clearScreen(RGB(250, 250, 250))
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN
				Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN
				sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), RGB(0, 0, 0))
				Dim As int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))
				Draw String (scrnPos.x + 2, scrnPos.y - 16), Str(sector(xi, yi).numInUse), 0
				'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0
			Next
		Next
		sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))
		For i As Integer = 0 To UBound(person)
			Dim As ULong c = stateColor(person(i).state)
			sg.drawCircleFilled(person(i).p, person(i).r, c)
		Next
		Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)
		Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)
		Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)
		Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)
		Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)
		For i As Integer = 1 To 9
			Draw String (10, 90 + i * 16), "Step: " & i & " = " &  Format(tAcc(i), "0.#0"), RGB(0, 0, 0)
		Next
		ScreenUnLock
		Sleep 1,1 'don't sleep every loop
	End If
	simulTime += timeStep
	simulDays = simulTime * DAY_PER_SEC
	tAcc(9) = 0 'sum other timers
	For i As Integer = 1 To 8
		tAcc(i) += (tNow(i) - tNow(i - 1))
		tAcc(9) += tAcc(i) 'sum
	Next
	If statCounters(ST_INFEC) = 0 Then Exit While
Wend
Locate 1,1: Print "End of simulation!" 
While InKey = "" : Sleep 1 : Wend
End

Next steps:
* Find optimal sector / grid size
* Make time graphs
* Apply social distancing measures

Edit: Current grid of 48 x 36 seems close to optimal. Larger grid is really slower, smaller grid does not give much improvement.

[badidea]

Simulation of virus containment failure:

Same conditions without borders:

Code: Select all

#Include "fbgfx.bi"
#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

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

'list can grow, but never shrink, for performance, non-sorted
Type grow_list
	Dim As Integer index(Any)
	Private:
	Dim As Integer numItems
	Public:
	Declare Constructor(startSize As Integer)
	Declare Constructor()
	Declare Destructor()
	Declare Sub Add(newIndex As Integer)
	Declare Sub empty()
	Declare Function numAlloc() As Integer
	Declare Function numInUse() As Integer
	Declare Sub show()
	'non-list methods
End Type

Constructor grow_list(startSize As Integer)
	If startSize > 0 Then
		ReDim index(startSize - 1)
	End If
End Constructor

Constructor grow_list()
	This.constructor(0)
End Constructor

Destructor grow_list()
	Erase(index)
End Destructor

Sub grow_list.add(newIndex As Integer)
	Dim As Integer ub = UBound(index)
	'if list is full, increase list size by 1
	If numItems = ub + 1 Then
		ReDim Preserve index(numItems)
	End If
	index(numItems) = newIndex
	numItems += 1
End Sub

Sub grow_list.empty()
	numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
	Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
	Return numItems
End Function

'for debugging
Sub grow_list.show()
	Print "--- " & numInUse() & " / " & numAlloc() & " ---"
	For i As Integer = 0 To numItems - 1
		Print i, index(i)
	Next
End Sub

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

#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

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

'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 drawCircleFilled(p As sgl2d, r As Single, c 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)
	Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)
	Declare Sub drawRectFilled(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)
	Dim As int2d posScrn = pos2screen(p)
	Circle(posScrn.x, posScrn.y), r * scale, c,,,,f
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

Sub scaled_graphics_type.drawRect(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, b
End Sub

Sub scaled_graphics_type.drawRectFilled(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, bf
End Sub

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

Const As Single INFEC_DIST = 0.05 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.25 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY 

Const PERSONS = 5000
Const As Single V_MAX = 0.1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10 
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)

'screen stuff
Const As Single PPM = 6.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.setScreen(SCRN_W, SCRN_H)
sg.setScaling(PPM, sgl2d(0, 0))

Enum HEALTH_STATE
	ST_INIT '0
	ST_INFEC '1
	ST_RECOV '2
	ST_DEAD '3
	ST_LAST 'number of states
End Enum

Dim As Integer statCounters(ST_LAST - 1)
Dim As ULong stateColor(ST_LAST - 1) = _
	{RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}
Dim As Single statPercentage(ST_LAST - 1)

Const As ULong WHITE = RGB(250, 250, 250)
Const As ULong BLACK = RGB(0, 0, 0)

Type person_type
	p As sgl2d 'position [m]
	pPrev As sgl2d 'position [m]
	v As sgl2d 'velocity [m/s]
	r As Single 'radius [m]
	sickEndTime As Single 'not double!
	state As Integer 'health
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 Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
	If value < min Then value = min
	If value > max Then value = max
End Sub

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

Randomize 1234
'initialise persons
For i As Integer = 0 To UBound(person)
	person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm
	person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm
	person(i).pPrev = person(i).p
	person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
	person(i).r = 0.35
	'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))
	person(i).state = ST_INIT
Next

'time step such that a max speed, position change = 20 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.2) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

'move view
sg.offset.x = -35
sg.offset.y = -25
Dim As Double tNow(0 To 9), tAcc(0 To 9)
Dim As Integer trigger = 1

sg.clearScreen(WHITE)
While Not MultiKey(FB.SC_ESCAPE)
	tNow(0) = Timer

	'trigger outbreak
	If MultiKey(FB.SC_SPACE) Then trigger = 1
	If trigger = 1 Then
		trigger = 0
		person(0).state = ST_INFEC
		person(0).sickEndTime = simulTime _
			+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
	End If
	'zoom / drag view by mouse
	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
	'patient 0 perpective
	If FOLLOW_PATIENT_0 = TRUE Then
		sg.offset.x = (person(0).p.x)
		sg.offset.y = (person(0).p.y)
	End If

	tNow(1) = Timer

	'update positions
	For i As Integer = 0 To UBound(person)
		person(i).p += person(i).v * timeStep
	Next
	'random direction change
	For i As Integer = 1 To 1
		Dim As Integer iPerson = Int(rndRangeSgl(0, PERSONS))
		If Rnd() > 0.5 Then
			person(iPerson).v.x = -person(iPerson).v.x
		Else
			person(iPerson).v.y = -person(iPerson).v.y
		End If
	Next

	tNow(2) = Timer

	'check wall collisions
	For i As Integer = 0 To UBound(person)
		With person(i)
			If .p.x < MAP_X_MIN Then .v.x = +Abs(.v.x)
			If .p.x > MAP_X_MAX Then .v.x = -abs(.v.x)
			If .p.y < MAP_Y_MIN Then .v.y = +Abs(.v.y)
			If .p.y > MAP_Y_MAX Then .v.y = -abs(.v.y)
			If .pPrev.x > 0 And .p.x < 0 Then .v.x = +Abs(.v.x)
			If .pPrev.x < 0 And .p.x > 0 Then .v.x = -abs(.v.x)
			If .pPrev.y > 0 And .p.y < 0 Then .v.y = +Abs(.v.y)
			If .pPrev.y < 0 And .p.y > 0 Then .v.y = -abs(.v.y)
		End With
	Next

	tNow(3) = Timer

	'check end of sickness
	For i As Integer = 0 To UBound(person)
		If person(i).state = ST_INFEC Then
			If simulTime > person(i).sickEndTime Then
				If Rnd() < MORTALITY Then
					person(i).state = ST_DEAD
					person(i).v = sgl2d(0, 0)
				Else
					person(i).state = ST_RECOV
					'person(i).v /= SICK_SPEED_FACTOR 'somthing wrong with this
					'person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
					'person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s
				End If
			End If
		End If
	Next
	
	tNow(4) = Timer

	'clear sectors
	For xi As Integer = 0 To NUM_SECT_X - 1
		For yi As Integer = 0 To NUM_SECT_Y - 1
			sector(xi, yi).empty()
		Next
	Next
	
	tNow(5) = Timer

	'assign persons to sectors
	For i As Integer = 0 To UBound(person)
		Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))
		Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))
		limitInt(xi, 0, NUM_SECT_X - 1)
		limitInt(yi, 0, NUM_SECT_Y - 1)
		sector(xi, yi).add(i)
	Next

	tNow(6) = Timer

	'spread the virus
	If sectioning = TRUE Then
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				'loop source in 1 sector
				For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1
					Dim As Integer iSrc = sector(xi, yi).index(iiSrc)
					If person(iSrc).state = ST_INFEC Then
						'check against targets in near sectors, including this sector
						For xdi As Integer = -1 To +1
							If xi + xdi < 0 Then Continue For
							If xi + xdi >= NUM_SECT_X Then Continue For
							For ydi As Integer = -1 To +1
								If yi + ydi < 0 Then Continue For
								If yi + ydi >= NUM_SECT_Y Then Continue For
								'loop targets with 1 (near) sector
								For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1
									Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)
									If person(iTar).state = ST_INIT Then
										If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
											person(iTar).state = ST_INFEC
											person(iTar).sickEndTime = simulTime _
												+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
											'person(iTar).v *= SICK_SPEED_FACTOR
										End If
									End If
								Next
							Next
						Next
					End If
				Next
			Next
		Next
	Else
		'case: sectioning = FALSE
		'check for disease transmission (iSrc = source, iTar = target)
		For iSrc As Integer = 0 To UBound(person)
			If person(iSrc).state = ST_INFEC Then
				For iTar As Integer = 0 To UBound(person)
					If person(iTar).state = ST_INIT Then
						If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then
							person(iTar).state = ST_INFEC
							person(iTar).sickEndTime = simulTime _
								+ rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY
							'person(iTar).v *= SICK_SPEED_FACTOR
						End If
					End If
				Next
			End If
		Next
	End If

	tNow(7) = Timer

	'clear stats
	For i As Integer = 0 To UBound(statCounters)
		statCounters(i) = 0
	Next
	'update stats
	For i As Integer = 0 To UBound(person)
		statCounters(person(i).state) += 1
	Next
	
	tNow(8) = Timer

	'draw world
	If (drawUpdate(DRAW_INTERVAL) = TRUE) Then
		ScreenLock
		'sg.clearScreen(WHITE)
		'clear part of screen
		sg.drawRectFilled(sgl2d(MAP_X_MIN - 2, MAP_Y_MIN - 2), _
			sgl2d(MAP_X_MAX + 2, MAP_Y_MAX + 2), WHITE)
		sg.drawLine(sgl2d(MAP_X_MIN, 0), sgl2d(MAP_X_MAX, 0), BLACK)
		sg.drawLine(sgl2d(0, MAP_Y_MIN), sgl2d(0, MAP_Y_MAX), BLACK)
		For xi As Integer = 0 To NUM_SECT_X - 1
			For yi As Integer = 0 To NUM_SECT_Y - 1
				Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN
				Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN
				Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN
				'sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), BLACK) 'grid
				'dim as int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))
				'draw string (scrnPos.x + 2, scrnPos.y - 16), str(sector(xi, yi).numInUse), 0
				'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0
			Next
		Next
		sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))
		For i As Integer = 0 To UBound(person)
			Dim As ULong c = stateColor(person(i).state)
			sg.drawCircleFilled(person(i).p, person(i).r, c)
		Next
		Line(0, 0)-(150, 100), WHITE, bf
		Draw String (10, 10 + 00), "Day: " & Format(simulDays, "0.#0"), BLACK
		Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)
		Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)
		Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)
		Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)
		'~ for i as integer = 1 to 9
			'~ draw string (10, 90 + i * 16), "Step: " & i & " = " &  format(tAcc(i), "0.#0"), BLACK
		'~ next
		For i As Integer = 0 To ST_LAST - 1
			statPercentage(i) = 100 * (statCounters(i) / PERSONS)
			Circle(10 + 4 * simulDays, (SCRN_H - 10) - 3 * statPercentage(i)), 1, stateColor(i),,,,f
		Next
		ScreenUnLock
		Sleep 1,1 'don't sleep every loop
	End If
	simulTime += timeStep
	simulDays = simulTime * DAY_PER_SEC
	tAcc(9) = 0 'sum other timers
	For i As Integer = 1 To 8
		tAcc(i) += (tNow(i) - tNow(i - 1))
		tAcc(9) += tAcc(i) 'sum
	Next
	If statCounters(ST_INFEC) = 0 Then Exit While
Wend
Draw String (10, SCRN_W - 26), "End of simulation", BLACK
While InKey = "" : Sleep 1 : Wend
End

[badidea]

A different version. No more fancy math, all integers and random walk. And no collision detection. I set a location on the map as infected when an infected person is at that spot. When a healthy (not immune) person is also at this spot, the disease is transmitted. Map is cleaned each cycle. Simulation is with 2500 persons, but 10 times more is no problem (switch display to pixel instead of circle). The daily movement of persons per day is very limited. I think that this shows that although 1 person is 'only' sick for 14 days, the total duration of the outbreak lasts much longer.

Code: Select all

const SCRN_W = 1000, SCRN_H = 700

screenres SCRN_W, SCRN_H, 32
width SCRN_W \ 18, SCRN_H \ 16

const ST_INIT = 0, ST_INFEC = 1, ST_RECOV = 2, ST_DEAD = 3, ST_INVALID = 4
const SICK_DAYS = 14
const as single MORTALITY = 3 / 100

const as ulong WHITE = rgb(200, 200, 200), BLACK = rgb(0, 0, 0)
dim as ulong stateColor(ST_INIT to ST_DEAD) = _
	{rgb(0, 200, 0), rgb(250, 0, 0), rgb(200, 200, 0), rgb(200, 200, 200)}
dim as integer statCounters(ST_INIT to ST_DEAD)
dim as string statNames(ST_INIT to ST_DEAD) = _
	{"Initial", "Infected", "Recovered", "Dead"}

const SEC_MIN = 60, SEC_HR = SEC_MIN * 60, SEC_DAY = SEC_HR * 24

type person_type
	as integer x, y
	as integer state
	as integer sickEndTime
end type

const PERSONS = 2500
dim as person_type person(any)
redim person(PERSONS - 1)

'map array dynamic, else too large for stack
const MAP_X = 700, MAP_Y = MAP_X
dim as integer map(any, any)
redim map(MAP_X - 1, MAP_Y - 1)

sub clearMap(map() as integer)
	clear map(0), 0, MAP_Y * MAP_X
	'set map borders
	for xi as integer = 0 to MAP_X - 1
		map(xi, 0) = ST_INVALID
		map(xi, MAP_Y - 1) = ST_INVALID
	next
	for yi as integer = 0 to MAP_Y - 1
		map(0, yi) = ST_INVALID
		map(MAP_X - 1, yi) = ST_INVALID
	next
end sub

clearMap(map())

const RND_DEF = 0, RND_CRT = 1, RND_FST = 2, RND_MTW = 3, RND_QBC = 4, RND_SYS = 5
randomize 222, RND_FST

for i as integer = 0 to PERSONS - 1
	person(i).x = int(rnd * (MAP_X - 2)) + 1
	person(i).y = int(rnd * (MAP_Y - 2)) + 1
next
person(0).state = ST_INFEC
person(0).sickEndTime = SICK_DAYS * SEC_DAY

dim as integer quit = 0
dim as integer simTime = 0, simTimeStep = SEC_MIN * 10
while quit = 0
	if inkey = chr(27) then quit = 1
	'move persons
	for i as integer = 0 to PERSONS - 1
		'dead don't walk
		if person(i).state < ST_DEAD then
			select case int(rnd * 4)
			case 0 'try right
				if map(person(i).x + 1, person(i).y) <> ST_INVALID then person(i).x += 1
			case 1 'try left
				if map(person(i).x - 1, person(i).y) <> ST_INVALID then person(i).x -= 1
			case 2 'try down
				if map(person(i).x, person(i).y + 1) <> ST_INVALID then person(i).y += 1
			case 3 'try up
				if map(person(i).x, person(i).y - 1) <> ST_INVALID then person(i).y -= 1
			end select
			'set infected locations
			if person(i).state = ST_INFEC then map(person(i).x, person(i).y) = ST_INFEC
		end if
	next
	'set new infections
	for i as integer = 0 to PERSONS - 1
		if person(i).state = ST_INIT then
			if map(person(i).x, person(i).y) = ST_INFEC then
				person(i).state = ST_INFEC
				person(i).sickEndTime = simTime + SICK_DAYS * SEC_DAY
			end if
		end if
	next
	'reset map & check end of sickness
	for i as integer = 0 to PERSONS - 1
		if person(i).state = ST_INFEC then
			map(person(i).x, person(i).y) = 0
			if simTime > person(i).sickEndTime then
				person(i).state = iif(rnd() < MORTALITY, ST_DEAD, ST_RECOV)
			end if
		end if
	next
	'clearMap(map())
	'clear stats
	for i as integer = 0 to ubound(statCounters)
		statCounters(i) = 0
	next
	'update stats
	for i as integer = 0 to ubound(person)
		statCounters(person(i).state) += 1
	next
	'draw
	screenlock
	line(0, 0) - (SCRN_W - 1, SCRN_H - 1), BLACK, bf
	for i as integer = 0 to PERSONS - 1
		'line(150 + person(i).x * 2, person(i).y * 2)-step(1, 1), stateColor(person(i).state), bf
		'pset(150 + person(i).x, person(i).y), stateColor(person(i).state)
		circle(150 + person(i).x, person(i).y), 3, stateColor(person(i).state),,,,f
	next
	locate 1, 1: print "days: " & (simTime \ SEC_DAY)
	for i as integer = 0 to ubound(statCounters)
		locate 2 + i, 1: print statNames(i) & ": " & statCounters(i)
	next
	screenunlock
	
	simTime += simTimeStep
	sleep 1
	if statCounters(ST_INFEC) = 0 then quit = 1
wend
sleep

[ShawnLG]

I have updated my CoronaSim with long distance vacation travel, simulate communal movement and use hashing to improve the simulator's time complexity from O(n^2) to O(n). Code is updated in the first post. It is like watching the Game of Life but with virus VS humans.

[badidea]

I have updated my CoronaSim with long distance vacation travel, simulate communal movement and use hashing to improve the simulator's time complexity from O(n^2) to O(n). Code is updated in the first post. It is like watching the Game of Life but with virus VS humans.

I have some weird lines flashing on the screen, what are those? Oh, that is the long distance travel.

[ShawnLG]

I have some weird lines flashing on the screen, what are those? Oh, that is the long distance travel.

Yes they are.