/* 
 *
 * Example-program: cml.c    1 dim coupled map lattice
 *
 */

# include <stdio.h>  

# include "nn.h"
# include "vector.h" /* matrix and vector math. */
# include "random.h" /* random numbers used     */


/* 
 *
 *  define simulation variables and parameters 
 *
 */



# define N 128      /* number of cells             */
Vector u;           /* state variables             */
Vector fu;          /* outputs f(x)                */
bVector y;          /* display variables (f(x)>1)  */
float m;            /* order parameter             */

/*
 *
 * define switch and slider variables needed
 * to build the grapical interface
 *
 */

SwitchValue ssynch; /* synch or asynch update */
SliderValue seps;   /* coupling strength      */
SliderValue sr;     /* parameter of f(x)      */



/* 
 *
 * Definition of graphical interface
 *
 */


BEGIN_DISPLAY

 SWITCH( "synch/asynch", ssynch) 
 SLIDER( "epsilon", seps, 0, 100)
 SLIDER( "r", sr, 0, 1000)

 WINDOW("dynamics")
  RASTER( "cml", AR, AC, y, bVECTOR, N, 0, .0, 1.0, 2)

 WINDOW("order param")
  GRAPH( "m(t)", AR, AC, &m, FLOAT_TYPE, 0, 0, 0, 0, -.05, 1.05 ) 

END_DISPLAY


NO_OUTPUT


float f( x )
float x;
{
  if (x > 0 && x < .5)
    return(.01 * sr * x);
  else if ( x >.5 && x < 1. )
    return(.01 * sr * (1.- x ) );
  else if (x > 1. && x < .005 * sr )
    return x;
  else
    return 0.; /* for appropriate initial conditions 
                               this should  not happen */ 
}

int main_init()
{
  u  = Get_Vector( N ); 
  fu = Get_Vector( N ); 
  y  = Get_bVector( N ); 
  randomize( 100 );
}

int init()
{
  int i;

  Clear_Vector( N, u);
  Clear_Vector( N, fu);
  Clear_bVector( N, y);
  for (i=0; i<N; i++)
  {
    u[i] = 2*equal_noise();
    y[i] = u[i] > 1.0 ? 0 : 1; 
  }
  m = bSum( N, y)/(float)N;
}

int step()
{
  int i,j;

  if (ssynch)  /* synchronous update */
  {
    for (j=0; j<N; j++)
      fu[j] = f(u[j]);
    for (j=0; j<N; j++)
    {
      u[j] = fu[j] + 0.005 * seps * ( fu[(j+N-1)%N] - 2*fu[j] + fu[(j+N+1)%N]);
      y[j] = u[j] > 1 ? 0 : 1;
    }
  }
  else         /* aynchronous update */
  {
    for (i=0; i<N; i++)
    {
      j = (int)(N*equal_noise());
      fu[j] = f(u[j]);
      u[j] = fu[j]+ 0.005 * seps * ( fu[(j+N-1)%N] - 2*fu[j] + fu[(j+N+1)%N]);
      y[j] = u[j] > 1 ? 0 : 1;
    }
  }
  m = bSum( N, y)/(float)N;
}