import array
import math
import random as r
import matplotlib.pyplot as plt
# from brachistochrone_fit import getCurve
# please email me for a program to fit the brachistochrone curve, the cycloid
# the program gives the theoretical curve, the cycloid for given endpoints

sqrt = math.sqrt

def dist(x1,y1,x2,y2): return sqrt((x1-x2)**2.0+(y1-y2)**2.0)
# the zero of energy is at h = 0, height h is measured downwards
def v(h,g): return sqrt(2.0*g*h) 
# the time required to get from 1 to 2
def Deltat(x1,y1,x2,y2,g):
    dh = y1 - y2
    if dh == 0.0:
        if v(y1,g) == 0.0: return float("infinity")
        return abs(x1-x2) / v(y1)
    else:
        return (v(y1,g)-v(y2,g)) * dist(x1,y1,x2,y2) / (g * dh)
def getLagrangian(X,Y,g):
    T = 0.0
    N = len(X)
    for i in range(1,N):
        T += Deltat(X[i],Y[i],X[i-1],Y[i-1],g)
    return T
def negative(y): return -y
    
# start coordinates
xs, ys = 0.0, 0.0
# end coordinates
xe, ye = 1.0, 0.0

# initial condition a direct slope
N = 50
X = array.array("f", [xs + i*(xe-xs)/float(N-1) for i in range(N)])
Y = array.array("f", [ys + i*(ye-ys)/float(N-1) for i in range(N)])

# gravity in natural units
g = 1.0
delta = 0.05
delta2 = 2 * delta
# inverse temperature, the fluctuations are smaller with lower temperatures
beta = 1500.0
IT= 2000000

print "starting run with: beta = ", str(beta)

for it in range(IT):
    # perturb a random particle, the ends are fixed
    n = r.randint(1,N-2)
    xp, yp = X[n] + delta2*r.random() - delta, Y[n] + delta2*r.random() - delta
    
    # to make sure there are no negative kinetic energies
    while (yp < 0.0):
        xp, yp = X[n] + delta2*r.random() - delta, Y[n] + delta2*r.random() - delta
    # change in travel time
    dt = Deltat(xp,yp,X[n-1],Y[n-1],g) \
        +Deltat(X[n+1],Y[n+1],xp,yp,g) \
        -Deltat(X[n],Y[n],X[n-1],Y[n-1],g) \
        -Deltat(X[n+1],Y[n+1],X[n],Y[n],g)
    P = math.exp(-beta*dt) if dt > 0 else 1.1
    p = r.random()
    if p < P: X[n], Y[n] = xp, yp
    if it % 1000000 == 0: 
        print "iteration: ", it, "L: ", getLagrangian(X,Y,g)

# please email me for a program to fit the brachistochrone curve, the cycloid
# the theoretical solution, the cycloid   
Xb = array.array('f', [0.0, 6.578438842552714e-06, 5.259635145193897e-05, 0.00017733754066284746, 0.00041977522778324783, 0.0008184179314412177, 0.0014111577766016126, 0.0022351210936903954, 0.0033265219535678625, 0.0047205169685184956, 0.006451072171330452, 0.008550820872187614, 0.011050945147871971, 0.013981038704514503, 0.01736900955438614, 0.021240947768092155, 0.025621041655540466, 0.03053145669400692, 0.035992302000522614, 0.04202147573232651, 0.04863465949892998, 0.055845193564891815, 0.06366412341594696, 0.07210004329681396, 0.08115911483764648, 0.09084506332874298, 0.10115910321474075, 0.1121000349521637, 0.12366412580013275, 0.13584521412849426, 0.148634672164917, 0.1620214730501175, 0.17599232494831085, 0.1905314326286316, 0.20562101900577545, 0.22124093770980835, 0.23736900091171265, 0.2539810240268707, 0.27105095982551575, 0.28855082392692566, 0.30645108222961426, 0.3247205317020416, 0.34332647919654846, 0.36223509907722473, 0.3814111351966858, 0.400818407535553, 0.4204197824001312, 0.4401773512363434, 0.46005260944366455, 0.480006605386734, 0.5, 0.5199934840202332, 0.5399473905563354, 0.5598226189613342, 0.5795801877975464, 0.599181592464447, 0.6185888648033142, 0.6377648711204529, 0.6566734910011292, 0.6752794981002808, 0.6935489773750305, 0.711449146270752, 0.7289490103721619, 0.7460189461708069, 0.7626310586929321, 0.7787591218948364, 0.794378936290741, 0.8094685077667236, 0.8240076303482056, 0.8379784822463989, 0.851365327835083, 0.8641547560691833, 0.876335859298706, 0.8878999352455139, 0.8988409042358398, 0.9091549515724182, 0.9188408851623535, 0.927899956703186, 0.9363358616828918, 0.9441547989845276, 0.9513653516769409, 0.9579785466194153, 0.9640077352523804, 0.9694685339927673, 0.9743789434432983, 0.9787590503692627, 0.9826309680938721, 0.98601895570755, 0.9889490604400635, 0.9914491772651672, 0.9935489296913147, 0.9952794909477234, 0.9966734647750854, 0.9977648854255676, 0.9985888600349426, 0.9991815686225891, 0.9995802044868469, 0.9998226761817932, 0.9999474287033081, 0.9999934434890747, 1.0])
Yb = array.array('f', [0.0, 0.00031405597110278904, 0.0012549844104796648, 0.0028190717566758394, 0.005000145640224218, 0.007789597846567631, 0.01117641944438219, 0.015147245489060879, 0.01968640461564064, 0.024775980040431023, 0.03039589151740074, 0.03652394935488701, 0.043135982006788254, 0.050205886363983154, 0.05770576372742653, 0.0656060203909874, 0.07387546449899673, 0.08248145878314972, 0.09139005839824677, 0.10056610405445099, 0.10997336357831955, 0.11957471072673798, 0.1293322741985321, 0.13920754194259644, 0.14916151762008667, 0.1591549515724182, 0.16914835572242737, 0.1791023463010788, 0.18897759914398193, 0.19873517751693726, 0.20833653211593628, 0.21774378418922424, 0.22691982984542847, 0.23582839965820312, 0.2444344162940979, 0.25270387530326843, 0.2606041133403778, 0.2681039869785309, 0.2751739025115967, 0.2817859351634979, 0.2879140079021454, 0.29353392124176025, 0.2986234724521637, 0.3031626343727112, 0.3071334660053253, 0.3105202913284302, 0.3133097290992737, 0.31549081206321716, 0.3170548975467682, 0.3179958164691925, 0.31830987334251404, 0.3179958164691925, 0.3170548975467682, 0.31549081206321716, 0.3133097290992737, 0.3105202913284302, 0.3071334660053253, 0.3031626343727112, 0.2986234724521637, 0.29353389143943787, 0.287913978099823, 0.2817859351634979, 0.2751739025115967, 0.26810401678085327, 0.2606040835380554, 0.25270384550094604, 0.24443446099758148, 0.2358284443616867, 0.22691984474658966, 0.21774379909038544, 0.20833654701709747, 0.19873517751693726, 0.18897761404514313, 0.1791023463010788, 0.16914837062358856, 0.15915493667125702, 0.14916151762008667, 0.13920752704143524, 0.1293322741985321, 0.11957470327615738, 0.10997334122657776, 0.1005660742521286, 0.09139003604650497, 0.08248143643140793, 0.07387549430131912, 0.06560604274272919, 0.05770578607916832, 0.05020590499043465, 0.04313599690794945, 0.0365239642560482, 0.03039589710533619, 0.02477598562836647, 0.01968640647828579, 0.015147245489060879, 0.011176418513059616, 0.007789595052599907, 0.005000142380595207, 0.002819068729877472, 0.0012549818493425846, 0.00031405442859977484, 2.4384253195449974e-15])

plt.title("g = %1.1f, T = %1.2f" % (g, getLagrangian(X,Y,g)))
plt.xlabel("x")
plt.ylabel("y")
plt.xlim([xs,xe])
plt.plot(X,map(negative,Y))
plt.plot(Xb,map(negative,Yb))
plt.savefig("brachistochrone"+"_g_%1.1f_b_%1.1f_%1.1f_%1.1f.png" % (g,beta,xe,ye))
plt.show()

print getLagrangian(X,Y,g), getLagrangian(Xb,Yb,g)