#! /usr/bin/env python

from numpy import *
import Gnuplot#, Gnuplot.funcutils

"""Try to solve a diffusion problem numerically"""


def demo():
    g = Gnuplot.Gnuplot(debug=1)

    n = 0
    n_0 = 1
    while n < n_0 :
       n+=1
       x = arange(0+(10.*n/(1.*n_0)),10+(10.*n/(1.*n_0)), dtype='float')
       y = sin(x)
       #d = Gnuplot.Data(x,y,title="let's see", with_="lines")
       g.xlabel("x-label")
       g.ylabel("y-label")
       g('set data style lines')  #this is how I give gnuplot arbitrary commands
       g.title("This should move")
       g.plot(x,y)
       g.hardcopy('gp_print.ps', enhanced=1, color=1) #produce hardcopy
       g.reset()
    raw_input('Please press the return key to continue...\n')
    

def bound(location,T):
    if location == float():
       return 1.-0.5*(sin(64*pi*T))
    else:
       return 0.5+0.3*(sin(32*pi*T))

def alpha_time_dep(T):
    return 1.-0.5*(sin(4*pi*T))

def march(L,n,lo,hi,dt,alpha,T,u):
    up = zeros(n+1)
#    u  = zeros(n+1)
    x = linspace(0,L,n+1) # grid points in x direction
    dx = float(L)/float(n)
    #alpha
    up[1:n] = u[1:n] + alpha*dt/dx/dx*(u[2:n+1] - 2*u[1:n] + u[0:n-1])
    up[0] = bound(lo,T)
    up[n] = bound(hi,T)
    u = up
    return up,x


# when executed, just run:
if __name__ == '__main__':
    import time
    L = 1.
    n = 40
    alpha = 5. #heat conduction coefficient
    dt = 0.00005  #time increment

    g = Gnuplot.Gnuplot()    
#    u = zeros(n+1)
    lo = 0.
    hi = 1.
    u = linspace(bound(lo,0.), bound(hi,0.),n+1)
    i = 0
    T = 0.
    t0 = time.time()
    while T < .1:
       #t0 = time.time()
       up,x=march(L,n,lo,hi, dt,alpha,T,u)
       #g('quit')
       g.reset
       #g('set data style lines')
       g('set yrange [0:1.5]')
       #g('set title "This is where I\'d like the time %d shown"')
       #title = 'set title ' + repr(T*1000)
       #g('bfit = gprintf("b = %s*10^%S",T*1000)')
       #g('set title bfit')
       #while time.time() < t0 + 0.001:
       #   a = 0
       g.plot(up)
       u = up
       T += dt
    t1 = time.time()
    input('Please press the return key to continue...\n')
    print('This script took {} seconds to run'.format(t1-t0))