31

#!/usr/bin/env python
#Copyright E.V.Tikhonov 2013
#Distributed under the terms of the GNU General Public License v2

import numpy as np
import math
from math import pi
from math import e
from matplotlib import pyplot as pl
from matplotlib import animation

g=9.81 #gravity constant
R=2    #vessel radius

class Solution(object):
    """
    Class which implements analytical solution of differential equation
    m\frac{dV}{dt}=\frac{4}{3}\pi r^3g(\rho_p-\rho_m)-6\pi\eta rV
    for particle with density rho_p moving through medium with density
    rho_m under gravity force.

    """
    def __init__(self,x,v,a,t,r=0.3,rho_m=1.261,rho_p=11.34,eta=14.9):
        """
                Parameters
        ----------
        x : array, float
            length-N array of spatial coordinate
        v : array, float
            length-N array of particle instantaneous velocity
        a : array, float
            length-N array of particle acceleration
        t : array, float
            length-N array of evenly distributed time coordinate
        r : float
            particle radius

        rho_p: float
            particle density, g*cm**-3
            aluminium:    2.70
            iron:        7.874
            copper:     8.96
            lead:       11.34
            gold:        11.6
            uranium:        19.1
        rho_m: float
            medium density, g*cm**-3
            mercury:        13.534
            glycerol:        1.261
            sulphuric acid:    1.8356
            ethylene glycol:    1.1132
            water:            0.9982
        eta: float
            medium viscosity, P
            glycerol:         14.9
            supluric acid:    0.242
            ethylene glycol:    0.161
            mercury:        1.526e-2
            water:            8.94e-3
        """
        self.x,self.v,self.a,self.t = map (np.asarray,(x,v,a,t))
        N = self.x.size
        N = self.v.size
        N = self.a.size
        N = self.t.size
        self.r=r
        assert self.x.shape == (N,)
        assert self.v.shape == (N,)
        assert self.a.shape == (N,)
        assert self.t.shape == (N,)
        self.N = len(t)
        dt = t[1]-t[0]
#       A = (4/3.)*pi*r**3*g*(rho_p-rho_m)
        m = rho_p*(4/3.)*pi*r**3
        def set_v_from_t(_t):
            result=(2/9.)*(r**2*g*(rho_p-rho_m)/eta)*(1-e**-((6*pi*eta*r*_t)/m))
            return result;
        def set_a_from_t(_t):
            result=(4/3.)*(pi*r**3/m)*g*(rho_p-rho_m)*e**-((6*pi*eta*r*_t)/m)
            return result;
        _x=0
        for i in range(N):
            if _x >= 60:
                x[i]=60
                v[i]=0
                a[i]=0
            else:
                v[i] = float(set_v_from_t(t[i]))
                x[i] = _x+v[i]*dt
                a[i] = set_a_from_t(t[i])
                _x=x[i]
                print t[i],x[i],v[i],a[i]

N = 2**13
dt = 0.01
x = 0.0 * np.arange(N)
v = 0.0 * np.arange(N)
a = 0.0 * np.arange(N)
t = dt * (np.arange(N))
#t_max = 120
x_max=60
#N_steps = 50
frames = N
S = Solution(x=x,v=v,a=a,t=t)

fig=pl.figure()
hlim=(-R,R)
xmin=0
#xmax=S.x[S.N-1]
xmax=60
ax=fig.add_subplot(111,xlim=(-R*5,R*5),ylim=(xmax+xmax/10,xmin-xmax/10))
#ax.set_xticks((-R,0,R))
#ax.set_yticks((0,5,10,15))
pl.axhspan(0,60,xmin=0.5-0.05*R,xmax=0.5+0.05*R,facecolor='#8888ff',alpha=0.5)
title = ax.set_title("")
part, = ax.plot([],[],'o',ms=45*S.r)
time_text = ax.text(0, 0.9, '', transform=ax.transAxes)
time_template = 'time = %.2fs'
def init():
    part.set_data([],[])
    time_text.set_text('')
    return part, time_text

def animate(i):
    part.set_data(0.0,S.x[i])
    time_text.set_text(time_template%(S.t[i]))
    return part,time_text


anim=animation.FuncAnimation(fig, animate, init_func=init,frames=frames, interval=10, blit=True)

pl.show()