Python 041 Winkel zwischen landender Kugel und Oberflächennormale

Das Programm gibt den Winkel an der zwischen dem Vektor einer auf den Boden fallenden Kugel und der Oberflächennormale liegt.

from visual import curve, scene, sphere, color, array
from math import *
import numpy as np
import time

scene.center=(10,0,5)

l = 2.5
#coordinate system
curve(pos=[(-l/2,0,0),(2*l,0,0)],radius=0.05,color=color.green)
curve(pos=[(0,-l/2,0),(0,2*l,0)],radius=0.05,color=color.blue)
curve(pos=[(0,0,-l/2),(0,0,2*l)],radius=0.05,color=color.red)

xwidth=22
zwidth=11
ebene = []
normal = []

#ground surface
for x in range(xwidth):
ebene.append([])
normal.append([])
for z in range(zwidth):
ebene[x].append(sin(x/4.0)+sin(z/4.0)) #FUNCTION
normal[x].append(0)

#draw ground surface
for x in range(xwidth-1):
for z in range(zwidth-1):
curve(pos=[(x,ebene[x][z],z),(x+1,ebene[x+1][z],z)],radius=.1,color=color.green)
curve(pos=[(x,ebene[x][z],z),(x,ebene[x][z+1],z+1)],radius=.1,color=color.green)

#ball
ball = sphere(pos=(10,5,5),radius=.3)
vektor = [-0,-0.02,-0]

#calculate surface normal
def surfacenormal(pointx,pointz):
pointy = sin(pointx/4.0)+sin(pointz/4.0) #FUNCTION
start = (pointx,pointy,pointz)
#derivative and cross product of fw and fu
fu = np.array([1,0.25*cos(pointx/4.0),0]) #FUNCTION
fw = np.array([0,0.25*cos(pointz/4.0),1]) #FUNCTION
fv = np.cross(fw,fu)
end = (fv[0]+pointx,fv[1]+pointy,fv[2]+pointz)
normal = curve(pos=[start, end])
#print(normal.pos)
return normal.pos

t=0
veca = [0,0,0]
vecc = [0,0,0]
#movement
while 1:
t += 1
ball.pos = ball.pos + vektor
if (ball.pos.y <= (sin(int(ball.pos.x/4.0))+sin(int(ball.pos.z/4.0)))): #FUNCTION #maybe round()
veca[0] = -vektor[0] #incoming vector
veca[1] = -vektor[1]
veca[2] = -vektor[2]
vecb = surfacenormal(int(ball.pos.x),int(ball.pos.z)) #surface normal
vecc[0] = vecb[1][0]-vecb[0][0]
vecc[1] = vecb[1][1]-vecb[0][1]
vecc[2] = vecb[1][2]-vecb[0][2]
zaehler = veca[0]*vecc[0]+veca[1]*vecc[1]+veca[2]*vecc[2]
nenner = (np.sqrt((veca[0])**2+(veca[1])**2+(veca[2])**2)*np.sqrt((vecc[0])**2+(vecc[1])**2+(vecc[2])**2))
cosalpha = zaehler/nenner
alpha = np.arccos(cosalpha)
winkel = np.degrees(alpha)
print(winkel)

time.sleep(.01)