############Script to produce this figure (python)##############
from math import *
import numpy as np
import matplotlib.pyplot as plt

tarr=np.linspace(-1.0,3*pi,400)
xarr=np.zeros(len(tarr))
yarr=np.zeros(len(tarr))
vxarr=np.zeros(len(tarr))
vyarr=np.zeros(len(tarr))
axarr=np.zeros(len(tarr))
ayarr=np.zeros(len(tarr))
apxarr=np.zeros(len(tarr))
apyarr=np.zeros(len(tarr))
aZpxarr=np.zeros(len(tarr))
aZpyarr=np.zeros(len(tarr))
aFZp=np.zeros(len(tarr))
aF=np.zeros(len(tarr))
i=0
R=2.0
r=R
omega=1.0
for t in tarr:
	xarr[i]=R*omega*t-r*sin(omega*t)
	yarr[i]=R-r*cos(omega*t)
	i+=1

vxarr=np.diff(xarr)
vyarr=np.diff(yarr)

axarr=np.diff(vxarr)
ayarr=np.diff(vyarr)

#Kruemmungsradius
#vsquared
for i in range(len(tarr)-2):
	v2=vxarr[i]*vxarr[i]+vyarr[i]*vyarr[i]
	av=vxarr[i]*axarr[i]+vyarr[i]*ayarr[i]
	apxarr[i]=av*vxarr[i]/v2
	apyarr[i]=av*vyarr[i]/v2
      	aZpxarr[i]=axarr[i]-apxarr[i]
      	aZpyarr[i]=ayarr[i]-apyarr[i]
	aF[i]=sqrt(axarr[i]*axarr[i]+ayarr[i]*ayarr[i])
	aFZp[i]=sqrt(aZpxarr[i]*aZpxarr[i]+aZpyarr[i]*aZpyarr[i])


plt.figure()

plt.subplot(211,aspect='equal',adjustable='box',title='Reference frame "Street": velocity and acceleration')
plt.plot(xarr,yarr,'k')
plt.hold(1)
plt.quiver(xarr[::20],yarr[::20],vxarr[::20],vyarr[::20],color='r',label='velocity')
plt.hold(1)
plt.quiver(xarr[::20],yarr[::20],axarr[::20],ayarr[::20],color='g',label='acceleration')
plt.ylim([0,max(yarr)])
plt.xlim([0,max(xarr)])

#plt.axes().set_aspect('equal')
#pl.show()
#plt.axes().set_aspect('equal')
plt.subplot(212,aspect='equal',adjustable='box',title='Reference frame "Street": Tangent and Centripetal Force')
#plt.figure()
plt.plot(xarr,yarr,'k')
plt.hold(1)
plt.quiver(xarr[::20],yarr[::20],apxarr[::20],apyarr[::20],color='b')
plt.quiver(xarr[::20],yarr[::20],aZpxarr[::20],aZpyarr[::20],color='y')
plt.ylim([0,max(yarr)])
plt.xlim([0,max(xarr)])
plt.savefig('centripetal_force_vs_external_force_vectors.svg')
#plt.axes().set_aspect('equal')

#plt.show()
plt.figure()
plt.title("External Force, tangent and centripetal components")
plt.plot(axarr,ayarr,'g', label='External Force')
plt.hold(1)
plt.plot(aZpxarr,aZpyarr,'y.', label='Centripetal Force')
plt.plot(apxarr,apyarr,'b.',label='Tangent Force')
plt.ylim([min(min(ayarr),min(aZpyarr),min(apyarr)),max(max(ayarr),max(aZpyarr),max(apyarr))*1.6])
plt.legend()
plt.axes().set_aspect('equal')
plt.savefig('centripetal_force_vs_external_force_circles.svg')

plt.figure()
plt.title('External Force and centripetal component vs. time')
plt.hold(1)
plt.plot(tarr,aFZp,'y',label='centripetal component')
plt.plot(tarr,aF,'g',label='External Force')
plt.legend(loc='lower left')
plt.xlim([0,tarr[len(tarr)-3]])
plt.savefig('centripetal_force_vs_external_force_time_series.svg')
#plt.show()