# Include everything from the visual module
from visual import *
from random import *
from visual.graph import *

dt = 0.01 #dt is delta t and is 0.1 second if we are using MKS

g = vector(0,9.8,0)
# Turn off autoscale
autoscale = 0

# Give the property range to the scene object
scene.range = (5, 5, 5)




# Defining mass_1 as a sphere with coordinates, radius and color
mass_1 = sphere(pos=(0, 0, 0), radius = 0.1, color = color.green)

# mass_1 has a mass of 5 kg
mass_1.mass = 5

# mass_1 has a velocity of 2 (unit?)
mass_1.velocity = vector(3.76, 0, 0)

#Definere netforce
NetForce = mass_1.mass*-g

#Definere acceleration
NetA = NetForce/mass_1.mass

# Wall definition
wall_1 = box(pos = (0, -1, 0), size = (10, 0.01, 10), color = color.white)


rlabel=label(pos=(0,-1.0,0), text='position', xoffset=0, yoffset=-12,
         height=15, border=10,box=0)


xdist = []
xdist.append(mass_1.x)





#loop for 5 seconds:
t = 0 #start at time 0

maxShots = 10

nShots = 0

graphwindow = gdisplay (xtitle = 'range x', ytitle = 'N', ymax = maxShots/4, xmax = 2)
distanceX = ghistogram (bins=arange(0.1,2,0.01),color=color.red,accumulate = 1, average=1)
graphwindow2 = gdisplay(xtitle = 'range x', ytitle = 'range z', xmax = 2, ymax = 2)
hits = gdots (gdisplay = graphwindow2, color = color.green)





while nShots < maxShots:
    dvx = gauss (3.76, 0.11)
    

    
    mass_1.pos=(0,0,0)
    mass_1.velocity = vector(dvx,0,0)
    while mass_1.pos.y>-1:
        rate(75) #slows animation by 1/arg seconds
        mass_1.velocity += NetA*dt
        mass_1.pos=mass_1.pos + mass_1.velocity*dt
        t+=dt
        rlabel.text="r = ("+str(mass_1.x)+") i + ("+str(mass_1.y)+") j"

        #Stopper bolden når den rammer "gulvet"
        if mass_1.y < wall_1.y:
            mass_1.velocity = vector(0,0,0)
    else: print nShots, mass_1.pos
    xdist.append(mass_1.x)

    distanceX.plot(data=xdist)
    hits.plot(pos=(mass_1.pos.x,mass_1.pos.z))
      
    nShots = nShots + 1

middelX=sum(xdist)/maxShots

print middelX

icount=1
devX=0
devX2=0
while icount < maxShots:
    devX += (xdist[icount]-middelX)
    devX2 += (xdist[icount]-middelX)**2
    icount +=1

devX = devX/maxShots
devX2 = devX2/maxShots

sigma = pow(devX2,0.5)
errMeanX = sigma/(pow(maxShots,0.5))

print devX
print devX2
print sigma
print errMeanX