#ifndef __SLIDERBLOCK_H__
#define __SLIDERBLOCK_H__

/******************************************************************************
* SLIDERBLOCK.H                                                               *
*                                                                             *
* The header file for SLIDERBLOCK.C                                           *
* Various helpful (and a few necessary) funtions for sliderblock analysis.    *
* Notation and motivation are taken from "Numerical Simulation of Earth-      *
* quake Sequences" by John Rundle and David Jackson (1977).                   *
*                                                                             *
*                                                                             *
* James Holliday                                                              *
* University of California - Davis                                            *
*                                                                             *
******************************************************************************/

/* Custom macros ----------------------------------------------------------- */
#define MIN(X,Y) ({typeof(X)x_=(X);typeof(Y)y_=(Y);(x_< y_)?x_:y_;})
#define MAX(X,Y) ({typeof(X)x_=(X);typeof(Y)y_=(Y);(x_< y_)?y_:x_;})

#define RAND(LOW,HIGH) (LOW+(HIGH-LOW)*((double)rand()/RAND_MAX))

#define VABS(V) (sqrt((V[0]*V[0])+(V[1]*V[1])+(V[2]*V[2])))

#define M_SQRT3  1.73205080757
#define true     1
#define false    0
#define INFINITY 1E50

/* Define different types of geometry -------------------------------------- */
typedef enum {

  CHAIN_LATTICE = 1,      /* Default for 1D */
  SQUARE_LATTICE,         /* Default for 2D */
  TRIANGULAR_LATTICE,
  HEX_LATTICE,
  CUBIC_LATTICE           /* Default for 3D */

} Geometry ;

/* Set Flags (powers of 2 for binaray combinations) ------------------------ */
typedef enum {

  MASSIVE_BLOCKS     =   1,
  PERIODIC_BC        =   2,
  ROUGH_FAULT        =   4,  /* Default */
  SMOOTH_FAULT       =   8,
  FLEXIBLE_FAULT     =  16,
  HOMOGENEOUS_FAULT  =  32,
  R3_INTERACTION     =  64,
  RANDOM_DROP        = 128

} Flags ;


/* Custom typedefs --------------------------------------------------------- */
typedef double Vector[3];
typedef unsigned short bool;


/* Slider Block struct ----------------------------------------------------- */
typedef struct {

  Vector Pos;  /* Global position vector (initialized so spacing = 1.0) */

  double M;    /* Mass of block */

  Vector U;    /* Displacement vector (x,y,z) */

  Vector Vel;  /* Velocity of the driver block */

  double Kl;   /* Leaf spring constant (connected to driving block) */

  unsigned Nk; /* Num of interacting blocks (blocks in interaction radius)   */
  void** pB;   /* Pointer to each block within the interaction radius        */
  double* Kc;  /* Vector of coil spring constants for each interacting block */

  double B;    /* Static friction force (time dependent failure strength) */
  double D;    /* Dynamical friction force */

  double Bu;   /* "Ultimate" breaking strangth (for anelasticity) */
  double Z;    /* Drop in time dependent friction force (for anelasticity) */

  Vector S;    /* Instantaneous stress */

} Block;


/* Parameter struct -------------------------------------------------------- */
typedef struct {

  unsigned Nx;   /* Number of blocks in X dimension */
  unsigned Ny;   /* Number of blocks in X dimension */
  unsigned Nz;   /* Number of blocks in X dimension */

  unsigned N;    /* Nx * Ny * Nz                    */

  double sigma;  /* Friction unit                   */
  double lambda; /* "Stiffness"                     */
  double tau1;   /* Time constant                   */
  double tau2;   /* Time constant                   */

  double R;      /* Radius of interaction (units of block separation) */

  Geometry G;    /* Geometry of the grid            */
  Flags    F;    /* Initialization flags            */

} Params;


/* Function prototypes ----------------------------------------------------- */
void setParams(int argc, char *argv[], Params *P);

void initialize(Block *list[], const Params *P);
void deallocate(Block  list[], const Params *P);

bool saveState(const Block list[], const Params *P, double T, char name[]);
bool loadState(Block *list[], Params *P, double *T, char name[]);
bool saveStress(const Block list[], const Params *P, double T, char name[]);

void randomizeStress(Block *list[], const Params *P);

double driveBlock(Block *block, const Params *P, double T);

double calcStress(Block *block);

int cascade(Block* list[], unsigned N, const Params *P);

double calcEnergy(Block list[], unsigned N);

double nextFailure(Block list[], const Params *P);

/*****************************************************************************/
#endif