x=np.linspace(1,40,10000)
y=3+.5*x
xini, yini=1,3.5
plt.plot(x,y,label='vector1')
am=np.stack((x,y))
am[:,0]
xy,length=am.shape
xatom,yatom=2,3

xm=am[0,-1]
ym=am[1,-1]
Vector1=[xatom-xini,yatom-yini,0]
Vector2=[xm-xini,ym-yini,0]
unitb=Vector2/np.linalg.norm(Vector2)
impactvector=Vector1-unitb*np.dot(Vector1, unitb)

plt.plot(x,y)
plt.scatter(xatom,yatom)
b=np.linalg.norm(impactvector)
plt.plot(x,yatom+(impactvector[1]/impactvector[0])*(x-xatom))
plt.xlim(0,4)
plt.ylim(2,4)