#-*- coding: iso-latin-1 -*-
#
# HardenbergDA_1.py: patrones de biomasa del modelo de Hardenberg
# 
# marzo 2009

import scipy
from LDynamicArray import *
from visual import *

DISPLAYINTERVAL=1
DerivZero = 1
Periodic  = 2
NullFront = 3

ModeFront = 2    ## defino condiciones de frontera (de las 3 de arriba)


class HardenbergDA:

   # la función __init__ inicializa datos, estructuras, etc. del objeto
   # que se crea usando f = Hardenberg(parámetros)
   # en este caso define los valores de todos los
   # parámetros del modelo de Hardenberg
   def __init__(self, N=100, dx=1.0, eps=0.2, gamma=0.8, beta=0.7):
      self.Dn = DynamicArray((N, N), zmin=0.0, zmax=.50)  ## DOS DynamicArray
      self.Dw = DynamicArray((N, N), zmin=0.0, zmax=.50)  ## Dn y Dw para n y w
      self.n = 0.3*scipy.rand(N,N)  
      self.w = 0.3*scipy.rand(N,N) 
      self.Dn.Frontier(self.n, ModeFront)
      self.Dn.Frontier(self.w, ModeFront)
      self.N = N
      self.dx = dx
      self.L = N*dx
      self.mu = 0.2    ### parámetros del problema (Hardenberg)
      self.gama = 1.6
      self.sigma = 1.6
      self.alfa = 3.
      self.beta = 3.
      self.delta = 100
      self.rho = 1.5
      self.p = 0.3 
      self.nu = 300.
      self.t = 0.    ### inicializa tiempo

   def rhs(self):    ## lado dercecho del sistema Reacción-Difusión (Hardenberg)
      self.dn_dt = self.gama*self.w/(1.+self.sigma*self.w)*self.n \
                   - self.n*self.n - self.mu*self.n + self.Dn.Lap_9(self.n)
      self.dw_dt = self.p - (1.-self.rho*self.n)*self.w - (self.w*self.w)*self.n \
                   + self.delta*self.Dw.Lap_9(self.w) \
                   - self.delta*self.beta*self.Dw.Lap_9(self.n)

   def step(self, dt):    ## el "paso": evaluar lado derecho y actualizar
      self.rhs()  # lado derecho del sistema de ecuaciones
      self.Dn.Frontier(self.dn_dt, ModeFront)  # condiciones de frontera
      self.Dw.Frontier(self.dw_dt, ModeFront)
      self.n += dt * self.dn_dt  ## método de Euler a cada ecuación dif.
      self.w += dt * self.dw_dt
      clip(self.n, 0.0, 1.0)  ## limita valores entre cero y uno
      clip(self.w, 0.0, 1.0)
      self.Dn.set_zmax(max(ravel(self.n)))  ## establece máximos
      self.Dw.set_zmax(max(ravel(self.w)))

   def run(self, T, dt):  ### esta es la función del ciclo principal
      count=0
      t0 = self.t
      while self.t < t0 + T:  ## ciclo de tiempo
         if count%DISPLAYINTERVAL==0:  ## despliegue de matrices(DynamicArray)
            self.Dn.setTitle('Biomasa. t = %f' % self.t)
            self.Dn.display(self.n)
            self.Dw.setTitle('Humedad. t = %f' % self.t)
            self.Dw.display(self.w)
         self.step(dt)
         self.t += dt
         count +=1


def demo():
   print "Patrones de Biomasa"
   print "  Hardenberg"
   f = HardenbergDA(N=64, dx=1.)   ## f es un objeto de la clase HardenbergDA
   f.run(2000., 0.001)             ## llamado a la función "run" del objeto f
   return f

### programa principal:
demo()