/* Program to compute the Fig 3 in 

Matsumoto & Chua & Tanaka:
Simplest Chaotic Non-autonomous Circuit
Physical Review A, Vol 30 No 2 August 1984 p.1155-1157

by numerical integration of the set of non-linear ordinary differential 
equations (ODE) describing the non-linear electronic circuit.

Author of the program: Pavel Pokorny
E-mail: Pavel.Pokorny@vscht.cz

This program uses the library of numerical functions Easynum,
available at
ftp://147.33.5.11/pub/xpplot/easynum.tar.Z
http://staff.vscht.cz/~pokornp/easynum/


To run it 
(1) save this code in a file circuit.c
(2) save easynum library in a file easynum.c
(3) compile it by
$ cc circuit.c -lm
(4) run it by
$ a.out

To see the results type
$ a.out | xpplot -

*/
# include "easynum.c"
/************************************/
/* rhs is the right-hand-side function in the ODE dx/dt = f(x,t) */

void rhs (time,vect,pars, func,g,J,dim)
/* input */
 real time;
 vector 
  vect, /* the state vector */
  pars, /* parameters */
/* output */
 func,  /* the vector of derivatives */
 g;     /* not used here, just for compatibility with stiff systems */
 matrix J;  /* dtto */
 int dim;   /* dimension, here 2, set in main */
{
 real 
  R  = 60,        /* Ohm */
  L  = 100E-6,    /* H   */
  C1 = 0.1E-6,    /* F   */
  C2 = 400E-12,   /* F   */
  f0 = 700E3,     /* Hz  */
  E0 = 0.1,       /* V   */
  a  = (C2-C1)/(2*C1*C2),
  b  = (C2+C1)/(2*C1*C2),
  omega = 2 * Pi * f0,
  q  = vect[0],
  i  = vect[1],
  Uc = a * Abs(q) + b * q + E0,
  Emax = pars[0],
  E = Emax * Sin (omega*time);

  func[0] = i;
  func[1] = (E - Uc - R * i) / L;
} /* rhs */
/************************************/
int main ()
{
 real t = 0,  /* time */
  Ts = 1.0 / 700E3; 
/* sampling period; 
   if equal to forcing period, then the Poincare section is computed */
 vector v,pars;
 int i,imax=10000, /* number of output points */
  dim=2;

 pars[0] = 2; /* forcing amplitude 2 V */

 for (i=0;i<dim;i++) v[i] = 0; 
  /* set initial condition 
   If 0 is a stationary point of your system
   (which is usually not the case in non-autonomous systems like this) 
   use some non-zero initial condition. */

 for (i=0;i<imax;i++) {
  integdt (&t,Ts,v,rhs,pars,dim);
/* integrate the system rhs for Ts seconds 
   both input and output is the state vector v */

  printf ("%G %G\n", v[0]*1E9 /* nC */ ,v[1]* 1E3  /* mA */);

/* print the results to stdout in 2 columns
   use redirection to a file
$ program > results

To split the results into 2 files, then use
$ col1 < results > q
$ col2 < results > i

where col1 is a one-line program 
awk '{print $1}'
similarly col2 is 
awk '{print $2}'

These programs (stored e.g. in /usr/local/bin) split 
data from a multi-column file to separate files.
To glue them together, use
$ paste q i > gi
*/

 };
 return (0);
} /* main */
/************************************/