CoCalc Public FilesDave and Henry Final project 213.ipynb
Authors: David Vanderweg, Henry White
Views : 126
Compute Environment: Ubuntu 20.04 (Default)
In [2]:
import numpy as np
from matplotlib import pyplot as plt

In [3]:
#Generating a 2 by 100 array of all peg locations

pegLocations=np.zeros([2,100])
for i in range(0,10):
for j in range(0,10):
#Saving x-locations of all pegs
if i%2==0:
pegLocations[0,j+i*10]=10*j
if i%2!=0:
pegLocations[0,j+i*10]=10*j+5
#Saving all y-locations of pegs
pegLocations[1,(i*10):(i+1)*10]=86.66-(8.66*(i+1))

In [4]:
#Defining gravitational constant
g=981 #cm/s

#Defining x and y acceleration expressions
#Note that if statements are used to account for direction of drag

def ax(Vx,C):
if Vx<=0:
ans=C*Vx**2
else:
ans = -C*Vx**2
return ans

def ay(Vy,C):
if Vy<=0:
ans=-g+C*(Vy)**2
else:
ans=-g-C*(Vy)**2
return ans

In [5]:
h=0.001 #initializing time step used throughout code

#Defining velocity function using forward Euler's method
def Vel(V0x,V0y,C):
vx=V0x+h*ax(V0x,C)
vy=V0y+h*ay(V0y,C)
return vx,vy

In [6]:
#defining a function to model the position of the ball using forward Euler's method

def Pos(x0,Vx,y0,Vy,h):
x=x0+h*Vx
y=y0+h*Vy
return x,y

In [7]:
#Defining a simple function to model bouncing off pegs. This function outputs velocities vx and vy

def bounce(x,y,pegx,pegy,vx,vy):
xdist=x-pegx
ydist=y-pegy
#Accounting for case where disk hits peg straight on causing the velocity in y direction to 'flip'
if xdist==0:
return vx,-vy

#Accounting for the case where the disk hits the peg at an angle
else:
theta=np.arctan(ydist/xdist)

#introducing randomness of up to 10% to angle that disk bounces off the peg
highThresh=theta+(theta*0.1)
lowThresh=theta-(theta*0.1)
thetaRand=float(np.random.uniform(lowThresh,highThresh))

#Applying vector reflection matrix
reflMat=np.array(([np.cos(2*thetaRand),np.sin(2*thetaRand)],[np.sin(thetaRand*2),-np.cos(2*thetaRand)]))
velMat=np.array(([vx],[vy]))
VxRef,VyRef=-np.matmul(reflMat,velMat)[0],-np.matmul(reflMat,velMat)[1]

#Introducing random reduction of velocity of up to 5% upon hitting peg
VxRef2=float(np.random.uniform(VxRef-(VxRef*0.05),VxRef))
VyRef2=float(np.random.uniform(VyRef-(VyRef*0.05),VyRef))

return VxRef2,VyRef2

In [13]:
#Applying simp 1/3 method of integration

def integrate(y, h):

ind = 1
tot = y[0] + y[-1]

for i in y[1:-1]:
if (ind % 2) == 0:
tot += 2 * i

else:
tot += 4 * i
ind += 1


In [14]:
def Plinko(xstart,drag,totaltime):
#Setting initial conditions
t=0
N=int(totaltime*1000)
C=drag
y0= np.random.uniform(89.8,90.2)
x0=np.random.uniform(xstart-0.2,xstart+0.2)
prevVx=np.random.uniform(0,0.2)
prevVy=np.random.uniform(0,0.2)

#Plotting the pegs in blue
plt.plot(pegLocations[0],pegLocations[1],'.')

#Initializing position and velocity arrays so that values can be plotted
xArray=np.zeros(N)
yArray=np.zeros(N)
vxArray=np.zeros(N)
vyArray=np.zeros(N)

#The dummy variable below is simply used to ensure bounce function is not triggered in consecutive time iterations as that causes the ball to move in physically impossible ways.
dumVar=0

#Initializing variables used to calculate total distance travelled by the disk
dist=0
totaldist=0

#Looping through all possible time values
for i in range(0,N):
velx,vely=Vel(prevVx,prevVy,C)
x1,y1=Pos(x0,velx,y0,vely,h)
dist=float(np.sqrt((x1-x0)**2+(y1-y0)**2))  #Change in position of the disk
totaldist=dist+totaldist
x0=np.copy(x1)
y0=np.copy(y1)

#Converting these values from 1x1 arrays to floats so functions work properly
x0f=float(x0)
y0f=float(y0)
velxf=float(velx)
velyf=float(vely)

#If statements are used to ensure that disk will go to opposite end of board if x=0 or x=100
if x0f > 100:
x0f = x0f-100
if x0f< 0:
x0f = x0f+100

#Plotting position of disk in red
plt.plot(x0f,y0f,'ro')

#storing position and velocity values in an array so they can be plotted later
xArray[i]=x0f
yArray[i]=y0f
vxArray[i]=velxf
vyArray[i]=velyf

#Stops iterating through time if disk reaches bottom of board
if y0f<=0:
t=0.001*i
numIterations=i
#determining which bin the disk landed in
if x0f<=5.0 or x0f>95.0:
Bin=int(0)
else:
Bin=(x0f+5)//10

break

#Looping through pegLocations array to determine if disk is bouncing off of a peg or falling at any given moment

for j in range(0,np.size(pegLocations[0,:])):
if (np.sqrt(((x0f-pegLocations[0,j])**2)+(y0f-pegLocations[1,j])**2)<=3.0) and dumVar==0:

#If above conditions are met, bounce function will be applied. The 3.0cm threshold is based off the sum of the radius of the peg and the disk (2.5cm+0.5cm).

posPegx=pegLocations[0,j]
posPegy=pegLocations[1,j]
newVel=bounce(x0f,y0f,posPegx,posPegy,velxf,velyf)

prevVx=np.copy(newVel[0])
prevVy=np.copy(newVel[1])

dumVar=1

break

#If rebounding function is not called, ball will just move according to velocity function
else:
prevVx=np.copy(velx)
prevVy=np.copy(vely)

#If the for loop reaches the end of the pegLocations array without triggering the bounce function, the disk experiences 'falling motion'. We set dumVar to 0 in this case because it is now physically possible for the disk to fall and hit another peg so the bounce function can be called.
if j==(np.size(pegLocations[0,:]))-1:
dumVar=0

#Truncating x,y,vx, and vy arrays and initializing time array
xArrayt=xArray[0:numIterations]
yArrayt=yArray[0:numIterations]
vxArrayt=vxArray[0:numIterations]
vyArrayt=vyArray[0:numIterations]
timeArray=np.zeros(numIterations)
for i in range(0,numIterations):
timeArray[i]=i*0.001

#Using fft to calculate the maximally occurring frequency
freq = np.fft.fft(vyArrayt[1:])
maxfreqoccur = np.argmax(freq[1:150])
freqarea = abs(integrate(freq[1:200],1))

#Plotting position vs t, velocity vs t, and frequency
plt.figure(2)
plt.plot(timeArray,xArrayt,'-',label="x(t)")
plt.plot(timeArray,yArrayt,'-',label="y(t)")
plt.xlabel("time [s]")
plt.ylabel("position [m]")
plt.legend()

plt.figure(3)
plt.plot(timeArray,vxArrayt,'-',label="vx(t)")
plt.plot(timeArray,vyArrayt,'-',label="vy(t)")
plt.xlabel("time [s]")
plt.ylabel("velocity [m/s]")
plt.legend()

plt.figure(4)
plt.plot(abs(freq[1:150]))
plt.xlabel("Frequency of change of direction [Hz]")
plt.ylabel("Number of Occurences")

return t,Bin,totaldist,maxfreqoccur,freqarea

In [12]:
#Copying above plinko function but removing all plotting parts as they're unnecessary for histogram generation

#Setting initial conditions
numIterations=0
Bin=0
t=0
N=int(totaltime*1000)
C=drag
y0= np.random.uniform(89.8,90.2)
x0=np.random.uniform(xstart-0.2,xstart+0.2)
prevVx=np.random.uniform(0,0.2)
prevVy=np.random.uniform(0,0.2)

#Initializing position and velocity arrays so that values        can be plotted
xArray=np.zeros(N)
yArray=np.zeros(N)
vxArray=np.zeros(N)
vyArray=np.zeros(N)

#The dummy variable below is simply used to ensure bounce function is not triggered in consecutive time iterations as that caused the ball to move in physically impossible ways.
dumVar=0
dist=0
totaldist=0

#Looping through all possible time values
for i in range(0,N):
velx,vely=Vel(prevVx,prevVy,C)
x1,y1=Pos(x0,velx,y0,vely,h)
dist=float(np.sqrt(abs(x1**2-x0**2)+abs(y1**2-y0**2)))
totaldist=dist+totaldist
x0=np.copy(x1)
y0=np.copy(y1)

#Converting these values from 1x1 arrays to floats so functions work properly
x0f=float(x0)
y0f=float(y0)
velxf=float(velx)
velyf=float(vely)
vyArray[i]=velyf

#If statements are used to ensure that disk will go to opposite end of board if x=0 or x=100
if x0f > 100:
x0f = x0f-100
if x0f< 0:
x0f = x0f+100

#Stops iterating through time if disk reaches bottom of board
if y0f<=0:
t=0.001*i
numIterations=i
#determining which bin the disk landed in
if x0f<=5.0 or x0f>95.0:
Bin=int(0)
else:
Bin=(x0f+5)//10

break

#Looping through pegLocations array to determine if disk is bouncing off of a peg or falling at any given moment
for j in range(0,np.size(pegLocations[0,:])):
if (np.sqrt(((x0f-pegLocations[0,j])**2)+(y0f-pegLocations[1,j])**2)<=3.0) and dumVar==0:  #If these conditions are met, bounce function will be applied. The 3.0cm threshold is based off the sum of the radius of the peg and the disk (2.5cm+0.5cm).

posPegx=pegLocations[0,j]
posPegy=pegLocations[1,j]
newVel=bounce(x0f,y0f,posPegx,posPegy,velxf,velyf)

prevVx=np.copy(newVel[0])
prevVy=np.copy(newVel[1])

dumVar=1

break

#If rebounding function is not called, ball will just move according to velocity function
else:
prevVx=np.copy(velx)
prevVy=np.copy(vely)

#If the for loop reaches the end of the pegLocations array without triggering the bounce function, the disk experiences 'falling motion'. We set dumVar to 0 in this case because it is now physically possible for the disk to fall and hit another peg so the bounce function can be called.
if j==(np.size(pegLocations[0,:]))-1:
dumVar=0

timeArray=np.zeros(numIterations)
vyArrayt=vyArray[0:numIterations]

for i in range(0,numIterations):
timeArray[i]=i*0.001

#max frequency and area
freq = np.fft.fft(vyArrayt[1:])
maxfreqoccur = np.argmax(freq[1:150])
freqarea = abs(integrate(freq,1))

return t,Bin,totaldist,maxfreqoccur, freqarea


# For a Total Drag Coefficient of 0.5

In [15]:
test,Bin,td,maxf,farea=Plinko(50,0.5,5.0)
print("It takes ",test, "seconds for the disk to move through the plinko board.")
print("The disk lands in bin",Bin)
print("The total distance travelled by disk is", td,"cm")
print("The Area Under the Frequency Curve is",farea,"s^-2" )
print("Maximum Frequency Occurance is at",maxf,"Hz")

It takes 3.753 seconds for the disk to move through the plinko board. The disk lands in bin 5.0 The total distance travelled by disk is 123.11054192049028 cm The Area Under the Frequency Curve is 31173.038126919964 s^-2 Maximum Frequency Occurance is at 10 Hz
In [ ]:
#generating histograms

#Note that code will take approximately 50 seconds to do 20 'runs' and 80 seconds for 30 runs

#Initializing the number of runs
num=20
tArr=np.zeros(num)
binArr=np.zeros(num)
tdArr=np.zeros(num)
maxfreqArray=np.zeros(num)
freqAreaArray=np.zeros(num)
for i in range(0,num):

#plotting histograms for all output values

plt.figure(1)
plt.hist(tArr,bins=5)
plt.xlabel("time [s]")
plt.ylabel("frequency of occurence")

plt.figure(2)
plt.hist(binArr,bins=10,range=(0,10))
plt.xlabel("Bin Number")
plt.ylabel("frequency of occurence")

plt.figure(3)
plt.hist(tdArr,bins=5)
plt.xlabel("Total displacement [m]")
plt.ylabel("frequency of occurence")

plt.figure(4)
plt.hist(maxfreqArray,bins=5)
plt.xlabel("Maximum Frequency [Hz]")
plt.ylabel("frequency of occurence")

plt.figure(5)
plt.hist(freqAreaArray,bins=5)
plt.xlabel("Integrated Area under all frequencies [s^-2]")
plt.ylabel("frequency of occurence")



# For a Total Drag Coefficient of 0.1

In [12]:
test,Bin,td,maxf,farea=Plinko(50,0.1,5.0)
print("It takes ",test, "seconds for the disk to move through the plinko board.")
print("The disk lands in bin",Bin)
print("The total distance travelled by disk is", td,"cm")
print("The Area Under the Frequency Curve is",farea,"s^-2" )
print("Maximum Frequency Occurance is at",maxf,"Hz")

It takes 1.744 seconds for the disk to move through the plinko board. The disk lands in bin 8.0 The total distance travelled by disk is 133.25542101575226 cm The Area Under the Frequency Curve is 74718.38236147874 s^-2 Maximum Frequency Occurance is at 0 Hz
In [13]:
#generating histograms

#Note that code will take approximately 30 seconds to do 20 'runs' and 50 seconds for 30 runs

#Initializing the number of runs
num=30
tArr=np.zeros(num)
binArr=np.zeros(num)
tdArr=np.zeros(num)
maxfreqArray=np.zeros(num)
freqAreaArray=np.zeros(num)
for i in range(0,num):

#plotting histograms for all output values

plt.figure(1)
plt.hist(tArr,bins=5)
plt.xlabel("time [s]")
plt.ylabel("frequency of occurence")

plt.figure(2)
plt.hist(binArr,bins=10,range=(0,10))
plt.xlabel("Bin Number")
plt.ylabel("frequency of occurence")

plt.figure(3)
plt.hist(tdArr,bins=5)
plt.xlabel("Total displacement [m]")
plt.ylabel("frequency of occurence")

plt.figure(4)
plt.hist(maxfreqArray,bins=5)
plt.xlabel("Maximum Frequency [Hz]")
plt.ylabel("frequency of occurence")

plt.figure(5)
plt.hist(freqAreaArray,bins=5)
plt.xlabel("Integrated Area under all frequencies [s^-2]")
plt.ylabel("frequency of occurence")


Text(0, 0.5, 'frequency of occurence')

# For a Total Drag Coefficient of 0.05

In [14]:
test,Bin,td,maxf,farea=Plinko(50,0.05,5.0)
print("It takes ",test, "seconds for the disk to move through the plinko board.")
print("The disk lands in bin",Bin)
print("The total distance travelled by disk is", td,"cm")
print("The Area Under the Frequency Curve is",farea,"s^-2" )
print("Maximum Frequency Occurance is at",maxf,"Hz")

It takes 1.997 seconds for the disk to move through the plinko board. The disk lands in bin 2.0 The total distance travelled by disk is 175.30274951349813 cm The Area Under the Frequency Curve is 210229.8617088856 s^-2 Maximum Frequency Occurance is at 7 Hz