#!/usr/bin/python
# Computing geodesics between points on a plane
# Written by Ganesh Swami - 12-JUL-2007
# http://ergodicity.iamganesh.com

import numpy
from numpy import *
from math import *
from pylab import *


def main(N=32,dt=0.0001,T=40000):
    X = linspace (0, 1, N)
    dX = 1.0/N

    y = 2.0 * numpy.sin(X*2*math.pi) + 2*X

    y_x = gradient (y,dX)
    v = y_x / numpy.sqrt (1 + y_x**2)
    dif = gradient(v,dX)

    # Plot solution
    plot (X,2*linspace(0,1,N),'g.-')

    for i in range(0,T):
        y_x = gradient (y,dX)
        v = y_x / numpy.sqrt (1 + y_x**2)

        dif_1 = gradient(v,dX)
        dif_1 [0] = 0.0
        dif_1 [-1] = 0.0

        y_1 = y + dt*dif_1

        e = linalg.norm(y_1 - y)
        d = linalg.norm(dif_1 - dif)
        f = linalg.norm(dif)

        y = y_1;
        dif = dif_1;

        if e < 1e-6 or f < 1e-6:
            print "Breaking due to tolerance limit"
            break

        if mod(i,T/50.0) == 0:
            plot (X,y, 'b')
            print ('%05d\t%3.4f\t%3.4f\t%3.4f' % (i,f, e, d))

    plot (X,y,'r.-')

    title ('Gradient Descent for solution of a line')
    savefig ('energy.png')
    show ()

if __name__=="__main__":
    main(N=128, T=40000)