# lorenz-circuit.R
#
# Simulate the Lorenz circuit described in [1].
#
# K. Myneni, 2010-02-03, Creative Consulting for Research & Education
#
# Notes:
#
# 1) This script is written for the free, multi-platform, R computing environment.
#    See [2] to obtain R, and [3] for notes on using R. Simply load this script
#    using,
#
#        source("lorenz-circuit.R")
#
#    from within R, and two graphical windows displaying figures similar to
#    fig. 2 and fig. 3, from ref. [1] will be computed and displayed.
#
# 2) The add-on package for R, odesolve, must be installed. See ref. [3] for
#    instructions. You may also refer to the documentation provided with the
#    odesolve package to understand how to setup and solve ordinary differential
#    equations.
#
# References:
#
# [1] N. J. Corron, 2010, A Simple Circuit Implementation of a Chaotic Lorenz System,
#     http://ccreweb.org/documents/physics/chaos/LorenzCircuit3.html 
#
# [2] http://www.r-project.org
#
# [3] K. Myneni, Notes on Using R, http://ccreweb.org/documents/programming/R-notes.html
#

## Load the add-on R package, odesolve, for solving ordinary differential
## equations.

require("odesolve")

## Define the differential equations to be solved (eqns. 4 from ref. [1])

derivs <- function( t, x, p ) {
	xdot = p[1]*(-x[1] + x[2])
	ydot = (p[2]/3)*(3 - x[3])*x[1] - x[2]
	zdot = -p[3]*x[3] + x[1]*x[2]
	return( list( c( xdot, ydot, zdot ) ) )
}

## Set parameters used in [1]
sigma = 10
R = 30
b = 8/3

params = c( sigma, R, b )

## Pick some initial values of the variables, near the attractor

x0 = c( -1, 1, 0 )

## Set up an array of time values for the solution

times <- seq(0, 1000, by = 0.01)

## Solve the differential equations, and create arrays, t, x, y, z 
## from the solution. Skip past the transient region, where the system 
## is moving from its initial point onto the attractor.

out <- lsoda( x0, times, derivs, params )
n = length(out[,1])
range = 5000:7000
t = out[range, 1]
x = out[range, 2]
y = out[range, 3]
z = out[range, 4]

## Make a plot to show the time series of x, y, and z, similar to 
## figure 2 in [1].

par(mfrow=c(3,1))
plot(t, x, type="l", col="blue")
plot(t, y, type="l", col="red")
plot(t, z, type="l", col="green")

## Open a new graphics output window and plot two views of the
## attractor, as shown in figure 3 of [1], using more points from
## the output to show the structure of the attractor more clearly.

x11()
range2 = 5000:10000
t = out[range2, 1]
x = out[range2, 2]
y = out[range2, 3]
z = out[range2, 4]
par(mfrow=c(1,2))
plot(x, y, type="l", col="blue")
plot(x, z, type="l", col="red")