/*
This file is part of EZRide.
Copyright 2007 Dwight Barkley, 
Copyright 2010 Vadim Biktashev & Andrew Foulkes.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published
by the Free Software Foundation, either version 3 of the License,
or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "ezride.h"
#ifndef _EZSTEP_
#define _EZSTEP_

/* ---------------------------------------------------------------------
 * We only use Barkley's and FHN's models in EZRide. However, you can
 * adjust the kinetics however you wish for other models.
 *
 * NOTE: We only use an explicit scheme for FHN kinetics. 
 * --------------------------------------------------------------------- */

#if FHN

  #define U_KINETICS(u,v) ( (u)+dt_o_eps*( (u)-( (u)*(u)*(u)/3 ) -(v) ) ) 
  #define V_KINETICS(u,v) ( (v)+dteps*( (u)+bet-(gam*(v)) ))

#else 

  #define U_THRESHOLD(v)  ( ( one_o_a*(v) ) + b_o_a )
  #define G(u,v)          ( (u)-(v) )

 #if EXPLICIT   
     /* Explicit u-kinetics */

  #define U_KINETICS(u,uth) ( (u)+dt_o_eps*(u)*(one-(u))*((u)-(uth)) )

 #else          
     /* Implicit u-kinetics */
     /* The original (Physica 49D) scheme can be obtained by defining
      * F2m and F2p to be one */

  #define F1m(u,uth) ( dt_o_eps*(one-(u))*((u)-(uth)) )
  #define F1p(u,uth) ( dt_o_eps*(u)*((u)-(uth)) )
  #define F2m(u,uth) ( one + dt_o_2eps*((uth)-(u)*(u)) )
  #define F2p(u,uth) ( one + dt_o_2eps*(two*(u) -(uth)-(u)*(u)) )

  #define U_KINETICS(u,uth) (                                      \
          ( (u) < (uth) ) ?                                        \
          (u)/(one-F1m(u,uth)*F2m(u,uth) ) :                       \
          ((u)+F1p(u,uth)*F2p(u,uth))/(one+F1p(u,uth)*F2p(u,uth)) )

  #endif

  #define V_KINETICS(u,v) ( (v)+dt*G(u,v) )

#endif


/* ---------------------------------------- 
 * Defining the diffusion macros:           
 * U_DIFFUSION, V_DIFFUSION, ZERO_USED_SUMS 
 * ---------------------------------------- */

#if V_DIFF_ON        
     /* v is diffusing */
  #define U_DIFFUSION(index_1)     (dt_o_wh2   * sigma_u[index_1]) 
  #define V_DIFFUSION(index_1)     (dtDv_o_wh2 * sigma_v[index_1]) 
  #define ZERO_USED_SUMS(index_1)  sigma_u[index_1] = sigma_v[index_1] = zero;
#else
     /* v is not diffusing */
  #define U_DIFFUSION(index_1)     (dt_o_wh2 * sigma_u[index_1]) 
  #define V_DIFFUSION(index_1)     zero
  #define ZERO_USED_SUMS(index_1)  sigma_u[index_1] = zero;
#endif

/* -------------------------------- 
 * Defining the spatial-sum macros: 
 * ADD_TO_U_SUM, ADD_TO_V_SUM     
 * -------------------------------- */

#if NINEPOINT
/* 9-point Laplacian formula */
  #define ADD_TO_U_SUM(index,index_2) \
   {Real stupid_cc = u[index]; \
     sigma_u[(index_2)]  -= twenty*stupid_cc; \
     sigma_u[(index_2)+J_INC]  += four*stupid_cc; \
     sigma_u[(index_2)-J_INC]  += four*stupid_cc; \
     sigma_u[(index_2)+I_INC]  += four*stupid_cc; \
     sigma_u[(index_2)-I_INC]  += four*stupid_cc; \
     sigma_u[(index_2)+I_INC+J_INC] += stupid_cc; \
     sigma_u[(index_2)+I_INC-J_INC] += stupid_cc; \
     sigma_u[(index_2)-I_INC+J_INC] += stupid_cc; \
     sigma_u[(index_2)-I_INC-J_INC] += stupid_cc; \
   }
#else
/* 5-point Laplacian formula */
  #define ADD_TO_U_SUM(index,index_2) \
   {Real stupid_cc = u[index]; \
      sigma_u[(index_2)]  -= four*stupid_cc; \
      sigma_u[(index_2)+J_INC] += stupid_cc; \
      sigma_u[(index_2)-J_INC] += stupid_cc; \
      sigma_u[(index_2)+I_INC] += stupid_cc; \
      sigma_u[(index_2)-I_INC] += stupid_cc; \
   }
#endif         

#if !V_DIFF_ON   
/* v is not diffusing */
  #define ADD_TO_V_SUM(index,index_2) 
#else            
/* v is diffusing */
#if NINEPOINT
/* 9-point Laplacian formula */
  #define ADD_TO_V_SUM(index,index_2) \
   {Real stupid_cc = v[index]; \
     sigma_v[(index_2)]  -= twenty*stupid_cc; \
     sigma_v[(index_2)+J_INC]  += four*stupid_cc; \
     sigma_v[(index_2)-J_INC]  += four*stupid_cc; \
     sigma_v[(index_2)+I_INC]  += four*stupid_cc; \
     sigma_v[(index_2)-I_INC]  += four*stupid_cc; \
     sigma_v[(index_2)+I_INC+J_INC] += stupid_cc; \
     sigma_v[(index_2)+I_INC-J_INC] += stupid_cc; \
     sigma_v[(index_2)-I_INC+J_INC] += stupid_cc; \
     sigma_v[(index_2)-I_INC-J_INC] += stupid_cc; \
   }
#else          
/* 5-point Laplacian formula */
  #define ADD_TO_V_SUM(index,index_2) \
   {Real stupid_cc = v[index]; \
      sigma_v[(index_2)]  -= four*stupid_cc; \
      sigma_v[(index_2)+J_INC] += stupid_cc; \
      sigma_v[(index_2)-J_INC] += stupid_cc; \
      sigma_v[(index_2)+I_INC] += stupid_cc; \
      sigma_v[(index_2)-I_INC] += stupid_cc; \
  }
#endif 
#endif
#endif /*  _EZSTEP_  */