### Post by David Drake on Mar 17, 2020 11:03:22 GMT -5

Just as an interesting exercise, I created the code below to simulate the effects of different factors on the spread of an infectious disease. This isn't scientific and it wasn't designed by an epidemiologist, but it is interesting to observe what happens as you vary population density, how contagious it is, fraction of the population that stays put, the length of time someone can spread the disease, and so forth. I think it helps emphasize the value of not moving around.

It isn't pretty or efficient, but it works. Enjoy!

It isn't pretty or efficient, but it works. Enjoy!

`population=100 'population`

seeds=5 'number of starting sick

wellCount=population-seeds

transferRate=0.50'probability of contracting on contact

recoCount=0'number of sims that were sick but recovered

fractionStationary=0.75

daysInfectious=21

prompt "Population: ";population

prompt "Seeds: ";seeds

prompt "Stationary fraction: ";fractionStationary

prompt "Transfer on contact fraction: ";transferRate

prompt "Days infectious: ";daysInfectious

sickCount=seeds

'nomainwin

UpperLeftX=1

UpperLeftY=1

WindowWidth=800

WindowHeight=600

graphicbox #1.g, 1, 1, 500, 500

graphicbox #1.t1, 1, 501, 800, 100

graphicbox #1.t2, 501, 1, 300, 500

open "Sim" for graphics_nf_nsb as #1

#1.g "fill white;flush"

#1.t1 "fill darkblue;flush"

#1.t2 "fill black;flush"

#1 "trapclose [q]"

dim sim(500,10)'500 sims, 10 parameters each

'1=x location

'2=y location

'3=x velocity

'4=y velocity

'5=wellness parameter (0=well, 1=sick, 2=recovered)

'6=time since initial infection

'set initial paramaters

'random location

for a =1 to population

sim(a,1)=int(rnd(1)*450)+20

sim(a,2)=int(rnd(1)*450)+20

sim(a,3)=(0.5-rnd(rnd(1)))*5

if sim(a,3)<1 and sim(a,3)>-1 then

if sim(a,3)<1 then sim(a,3)=sim(a,3)-2

if sim(a,3)>1 then sim(a,3)=sim(a,3)+2

end if

sim(a,4)=(0.5-rnd(rnd(1)))*5

if sim(a,4)<1 and sim(a,4)>-1 then

if sim(a,4)<1 then sim(a,4)=sim(a,4)-2

if sim(a,4)>1 then sim(a,4)=sim(a,4)+2

end if

if a/population<fractionStationary then sim(a,3)=0:sim(a,4)=0

next a

'make sprites

#1.t2 "down;place 1 1;color white;backcolor white;boxfilled 10 10;place 5 5;backcolor black;color black;circlefilled 5;place 5 15;color green;backcolor green;circlefilled 5"

#1.t2 "getbmp well 1 1 8 20"

#1.t2 "color red;backcolor red;circlefilled 4"

#1.t2 "getbmp sick 1 1 8 20"

#1.t2 "color blue;backcolor blue;circlefilled 4"

#1.t2 "getbmp reco 1 1 8 20"

#1.g "getbmp bg 1 1 500 500"

#1.g "background bg"

#1.g "drawbmp well 10 10;drawbmp sick 40 40;drawbmp reco 80 80"

for a = 1 to population

#1.g "addsprite sim";a;" well sick reco"

#1.g "spritexy sim";a;" ";sim(a,1);" ";sim(a,2)

#1.g "spritemovexy sim";a;" ";sim(a,3);" ";sim(a,4)

next a

for a = 1 to seeds ' seed sick people with random times since infaction

sim(a,5)=1

sim(a,6)=rnd(1)*14*24

#1.g "spriteimage sim";a;" sick"

next a

#1.g "drawsprites"

timer 25, [animate]

wait

[animate]

scan

if sickCount<=0 then

waitCount=waitCount+1

if waitCount>50 then

timer 0

goto [endSim]

end if

end if

hours=hours+1

#1.g "drawsprites"

for a = 1 to population

if sim(a,6)>0 then 'decrement time sick

sim(a,6)=sim(a,6)-1

if sim(a,6)<=0 then

sim(a,5)=2 'move to reco

#1.g "spriteimage sim";a;" reco"

recoCount=recoCount+1

sickCount=sickCount-1

end if

end if

#1.g "spritecollides sim";a;" col$"

#1.g "spritexy? sim";a;" x y"

if col$<>"" then

if rnd(1)<0.5 then sim(a,3)=sim(a,3)*(-1) 'change direction

if rnd(1)<0.5 then sim(a,4)=sim(a,4)*(-1) 'change direction

col=val(left$(after$(col$,"sim"),2)) 'get the colliding sim

if sim(a,5)=1 and sim(col,5)=0 and rnd(1)<transferRate then 'move from well to sick

sim(col,5)=1:sim(col,6)=daysInfectious*24

#1.g "spriteimage sim";col;" sick"

wellCount=wellCount-1

sickCount=sickCount+1

else

if sim(a,5)=0 and sim(col,5)=1 and rnd(1)<transferRate then

sim(a,5)=1:sim(a,6)=daysInfectious*24 'move from well to sick

#1.g "spriteimage sim";a;" sick"

wellCount=wellCount-1

sickCount=sickCount+1

end if

end if

end if

if x>=480 or x<=10 then sim(a,3)=sim(a,3)*(-1) 'change direction at perimeter

if y>=480 or y<=10 then sim(a,4)=sim(a,4)*(-1) 'change direction at perimeter

#1.g "spritemovexy sim";a;" ";sim(a,3);" ";sim(a,4)

next a

if sickCount>maxSickCount then maxSickCount=sickCount

#1.t1, "place 10 40;\Time: ";hours;" Well: ";wellCount;" Sick: ";sickCount;" Reco: ";recoCount;" "

#1.t2 "discard;cls"

#1.t2, "place 10 100;backcolor white;color black;boxfilled 40 0"

#1.t2, "place 40 100;backcolor white;color black;boxfilled 70 0"

#1.t2, "place 70 100;backcolor white;color black;boxfilled 100 0"

#1.t2, "place 40 100;backcolor pink;color black;boxfilled 70 ";100-(100*maxSickCount/population)

#1.t2, "place 10 100;backcolor green;color black;boxfilled 40 ";100-(100*wellCount/population)

#1.t2, "place 40 100;backcolor red;color black;boxfilled 70 ";100-(100*sickCount/population)

#1.t2, "place 70 100;backcolor blue;color black;boxfilled 100 ";100-(100*recoCount/population)

wait

[endSim]

timer 0

print "Simulation complete"

print "Population: ";population

print "Seed sick: ";seeds

print "Stationary: ";fractionStationary

print "Tran rate: ";transferRate

print "Days infec: ";daysInfectious

print

print "Hrs of end: ";hours

print "% infected: ";int(recoCount/population*100);"%"

print "Max % sick: ";int(maxSickCount/population*100);"%"

print "------------------------"

wait

[q]

timer 0

unloadbmp "well"

unloadbmp "sick"

unloadbmp "reco"

#1.g "cls;discard"

#1.t1 "cls;discard"

#1.t2 "cls;discard"

close #1

end