/*  rros.c
 *
 *  either chaotic Roessler oscillators or
 *  undamped harm. osc. with reset;
 *  various coupling structures
 *
 */

# include <stdio.h>

# include <nn.h>     /* always include this     */
# include <vector.h> /* matrix and vector math. */
# include <random.h> /* random numbers are used */
# include <numerics.h>  /* Runge-Kutta */

# define ROSSLER 1  // set this to 1 if you wand Roessler oscillators
                    // otherwise undamped harmonic oscill. with reset are used


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


# define STEPSIZE .01
float t;

# define N 64       /* number of units     */
# define n  3       /* order of diff.system */

Vector domega; 

Vector  x;      /* x1 ... xN, y1 .... yN, z1 .... zN */
Vector  dxdt;     /* derivatives */


Matrix  J;         /* connections  (if not meanfield couplings)  */
                   /* diffusive or random ....                   */

Vector  cfields;   /* coupling fields; either meanfield or diffusive 
                      or random connectivity */
float xx1, yy1;

SwitchValue smean = ON;       /* mean field coupling */
SwitchValue sdiffusive = OFF; /* diffusive coupling */
SwitchValue swrand = OFF;    /* random connections */

SliderValue somega = 1000;
SliderValue sdelomega = 100;
SliderValue sepsilon = 100;
SliderValue sa = 150;


BEGIN_DISPLAY

 SWITCH( "mean", smean )
 SWITCH( "diffusive", sdiffusive )
 SWITCH( "random", swrand )

 SLIDER( "mean omega", somega, 500, 1500)
 SLIDER( "delta omega", sdelomega, 0, 500)
 SLIDER( "coupling strength", sepsilon, 0, 500)
 SLIDER( "a", sa, 0, 500)

 WINDOW("signals")

  RASTER( "x", AR, AC, x, VECTOR,  N, 0, 0.0, 1.0, 2) 

 WINDOW("MF-xy-plot")  /* yy1 DOESNT WORK ?????????? */

  PLOT("x-y", AR, AC, &xx1, VECTOR, 1, 0, 0, 0, -20., 20.,
                      &yy1, VECTOR, 1, 0, 0, 0, -20., 20.,  2 );

 WINDOW("xy-plot")

  PLOT("x-y", AR, AC, x, VECTOR, N, n, 0, 0, -20., 20.,
                      x, VECTOR, N, n, 0, 1, -20., 20.,  2 );

 WINDOW("x(t)")

  GRAPH( "x1", AR, AC, x, VECTOR, N, 0, 0, 0, -20, 20 )
  GRAPH( "x2", AR, NC, x, VECTOR, N, 0, 1, 0, -20, 20 )

  GRAPH( "y1", NR, C0, &x[N], VECTOR, N, 0, 0, 0, -20, 20 )
  GRAPH( "y2", AR, NC, &x[N], VECTOR, N, 0, 1, 0, -20, 20 )

  GRAPH( "z1", NR, C0, &x[2*N], VECTOR, N, 0, 0, 0, 0., 20 )
  GRAPH( "z2", AR, NC, &x[2*N], VECTOR, N, 0, 1, 0, 0., 20 )


 WINDOW("MF")

  GRAPH( "x1", AR, AC, &xx1, VECTOR, 1, 0, 0, 0, -20., 20.)
  GRAPH( "x2", AR, NC, &yy1, VECTOR, 1, 0, 0, 0, -20., 20.)

END_DISPLAY



NO_OUTPUT



int main_init()
{

  SET_STEPSIZE( STEPSIZE )
  randomize( time(NULL) );

  J   = Get_Matrix(N,N);
  x   = Get_Vector(N*n);
  dxdt= Get_Vector(N*n);
  domega = Get_Vector(N);
  cfields = Get_Vector(N);
}




int init()
{
  int i;

  Clear_Vector(N, domega);
  Clear_Vector(N, cfields);
  Clear_Vector(N*n, x);
  Clear_Vector(N*n, dxdt);
  t = 0.0;

  for (i=0 ; i<N; i++)
  {
    domega[i] = -.5+(1.*i)/N;
    cfields[i] = 0.0;
    x[i] = 4.;
    x[N+i] = 4;
  }

  Clear_Matrix( N,N, J);
  Make_Matrix( N, N, J, 1, 1);
}




void derivs(x,y,dfdx)
float x;
float *y;
float *dfdx;
{
  int i,j;
  float omeg,a;

  a = .001*sa;

  for (i=0;i<N;i++)
  {  
    if (ROSSLER)  /* original roessler */
    { 
      omeg = .001*(somega + sdelomega*domega[i]);
      dfdx[i]     = -omeg*y[N+i] - y[2*N+i] + cfields[i];
      dfdx[N+i]   = omeg*y[i] + a*y[N+i];
      dfdx[2*N+i] = .4+y[2*N+i]*(y[i]-8.5);
    } 
    else  /* antisymm. undamped harm.osc. */
    {
      omeg = .001*(somega + sdelomega*domega[i]);
      dfdx[i]     =    a*y[i] -omeg*y[N+i] - y[2*N+i] + cfields[i];
      dfdx[N+i]   =  omeg*y[i] +   a*y[N+i];
      dfdx[2*N+i] = .4+y[2*N+i]*(y[i]-8.5); 
    }
  }
}



int step()
{
  int i;
  float mf,epsfac;
  static float tlast=-1,phi1;


  rk4(x, dxdt, N*n, t, step_size, x, derivs);


  epsfac = .001*sepsilon;

  if(swrand)  /* different amplitude scaling in alternatives ... */
  {
    Mult( N, N, J, x, cfields );   /* good luck! first components of
                                      systems are first N vals of x */ 
    epsfac /= N;
  }
  else if(sdiffusive) /* open boundaries */
  {
    cfields[0] = x[1]-x[0]; 
    for(i=1;i<N-1;i++)
     cfields[i] = x[i+1]+x[i-1]-2*x[i];
    cfields[N-1] = x[N-2]-x[N-1];
  }
  else if (smean) /* mean field */
  {
    mf = Sum(N, x)/(float)N;
    for(i=0;i<N;i++)
     cfields[i] = mf;
  }
  else
  {
    for(i=0;i<N;i++)
     cfields[i] = 0.;
  }

  for(i=0;i<N;i++)   /* scale with coupling strength */
    cfields[i] *= epsfac;

  xx1 = Sum( N, x)/N;
  yy1 = Sum( N, &x[N])/N;

  t+=step_size;

}