from numpy import *
from matplotlib.pyplot import *
import sys as sys




beta =1/3.
V    =650.
k0   =325

ratio=array([0.093,0.14,0.4,0.6,0.8])
alpha=ratio**(V/k0)
gamma=3/(1-beta)
E0=V*(1-beta)**2/(6*k0*beta)
#f_1_10=ratio**(2*beta/(1-beta)+V/(k0*velocity))*velocity**(-3/(1-beta))*exp(-6*k0*beta*velocity/(V*(1-beta)**2))
peak=(-gamma+sqrt(gamma**2-4*log(alpha)/E0))/(2./E0)

print ratio, peak

figure()
plot(ratio,peak)
show()