Category : C Source Code
Archive   : DYNAM_S1.ZIP
Filename : YMAPS.C

/********************************** YMAPS ************************************/
/********************* (C) 1986,7,8 by JAMES A. YORKE ************************/




#include "yinclud.h"



/**************** MAPS and INITIALIZATIONS OF MAPS *******************/

mapCubic() { /* Cubic polynomial C Cubic in complex
plane rho*Z^3 + (C1+iC2)*Z^2 + (C3+iC4)*Z +
C5+iC6 */
int i;
double y_temp[6];
double X,
XX,
Y,
YY,
XXmYY,
XY,
XY2,
Fx,
Fy;
int base;

X = y[0];
Y = y[1];
XX = X * X;
YY = Y * Y;
XXmYY = XX - YY;
XY = X * Y;
XY2 = 2 * XY;
/* y_temp[0] = real part of polynomial: */
/* y_temp[0] = imaginary part of polynomial: */
/***** y_temp[0] = rho*X*(XX-3*YY)
+C1*XXmYY -C2*XY2
+C3*X - C4*Y
+C5;
y_temp[1] = rho*Y*(3*XX -YY)
+C2*XXmYY +C1*XY2
+C4*X + C3*Y
+C6; ************/

y_temp[0] = C5;
y_temp[1] = C6;
if(C4 != 0) {
y_temp[0] -= C4 * Y;
y_temp[1] += C4 * X;
}
if(C3 != 0) {
y_temp[0] += C3 * X;
y_temp[1] += C3 * Y;
}
if(C2 != 0) {
y_temp[0] -= C2 * XY2;
y_temp[1] += C2 * XXmYY;
}
if(C1 != 0) {
y_temp[0] += C1 * XXmYY;
y_temp[1] += C1 * XY2;
}
if(rho != 0) {
y_temp[0] += rho * X * (XX - 3 * YY);
y_temp[1] += rho * Y * (3 * XX - YY);
}


if(XX + YY >= 10000) { /* to prevent overflows */
y[0] = 100;
y[1] = 100;
dot = dots;
printf("trajectory is going to infinity ");
beep();
}
for(i = 0; i < num_lyap; i++) {
/* this is used both in Lyapunov number
computation and for the Newton method and
orbit continuation routines */
base = lyapzero + vec_dim * i;
/* Cauchy-Riemann equations: Fx = Gy, Fy = -Gx */
Fx = rho * 3 * (XXmYY)
+ 2 * (C1 * X - C2 * Y)
+ C3;
Fy = -rho * 6 * XY
- 2 * (C1 * Y + C2 * X)
- C4;
y_temp[base]
= Fx * y[base]/* partial of real part w.r.t. X */
+Fy * y[base + 1];/* partial of real part wrt Y */
y_temp[base + 1]
= -Fy * y[base]/* partial of imag part wrt X */
+Fx * y[base + 1];
}
for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
y[i] = y_temp[i];
return;
}

initcubic() {
vec_dim = 2; /* the dimension of the Lyapunov vectors =
phase space dim */
num_lyap = 0; /* the number of Lyapunov numbers to be
computed <= vec_dim */
lyapzero = 2; /* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
X_upper = 1.; /* x scale */
X_lower = -1.2;
Y_upper = 1.1; /* y scale */
Y_lower = -1.1;
y[eqn0 + 0] = 0.; /* initialize x coord of y1 */
y[eqn0 + 1] = 0.; /* initialize y coord of y1 */
setequal(0, 1, eqn1); /* reinitializes y[] = y1 */
/* Now initialize all constants that are used in the map even if the initial
value is 0 */
rho = 1.;
C1 = 0.;
C2 = 0.; /* even constants that are zero must be
initialized here since otherwise they will
be initialized = -9999. and the menu knows
that constants so initialized should not be
accessible from the menu */
C3 = 0.0;
C4 = 0.;
C5 =.38;
C6 = 0.;
map = mapCubic; /* note mapCubic is defined in this file prior
to its mention here so it does not have to
be declared in initCubic */
return;
}


rotor() { /* x->x+y then y->Cy+rho*sin(x) */
double moduloAB();
double s,
c;
int i;
double y_temp[6];
int base;

y_temp[0] = y[0] + y[1];
s = sin(y_temp[0]);
y_temp[1] = C1 * y[1] + rho * s;

if(num_lyap > 0)
c = cos(y_temp[0]);
for(i = 0; i < num_lyap; i++) {
base = lyapzero + vec_dim * i;
y_temp[base]
= y[base]
+ y[base + 1];
y_temp[base + 1]
= rho * c * (y[base] + y[base + 1])
+ C1 * y[base + 1];
}
for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
y[i] = y_temp[i];

y[0] = moduloAB(y[0], -pi, pi);
if(C1 == 1) /* area preserving case */
y[1] = moduloAB(y[1], -pi, pi);
}


mapComplex() { /* Tinkerbell -- this map is not analytic */
int i;
double y_temp[6];
double X,
XX,
Y,
YY;
int base;
X = y[0];
Y = y[1];
XX = X * X;
YY = Y * Y;
y_temp[0] = XX - YY + C1 * X + C2 * Y;
y_temp[1] = 2 * X * Y + C3 * X + C4 * Y;
/* if(XX + YY >= 100000000) /to prevent overflows /
{y[0] = 100;
y[1] =100;
dot = dots;
printf("trajectory is going to infinity ");
beep();
}
*/
for(i = 0; i < num_lyap; i++) {
/* this is used both in Lyapunov number
computation and for the Newton method and
orbit continuation routines */
base = lyapzero + vec_dim * i;
y_temp[base]
= (2 * X + C1) * y[base]
+ (-2 * Y + C2) * y[base + 1];
y_temp[base + 1]
= (2 * Y + C3) * y[base]
+ (2 * X + C4) * y[base + 1];
}
for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
y[i] = y_temp[i];
return;
}





initComplex() {
vec_dim = 2; /* the dimension of the Lyapunov vectors =
phase space dim */
num_lyap = 0; /* the number of Lyapunov numbers to be
computed <= vec_dim */
lyapzero = 2; /* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
X_upper =.5; /* x scale */
X_lower = -1.3;
Y_upper = -1.6; /* y scale */
Y_lower = 0.6;
y[eqn0 + 0] = -.72;
y[eqn0 + 1] = -.64;
setequal(0, 1, eqn1); /* reinitializes y[] */
C1 =.9;
C2 = -.6013;
C3 = 2.0;
C4 =.5;
map = mapComplex;
}
initRotor() {
vec_dim = 2; /* the dimension of the Lyapunov vectors =
phase space dim */
num_lyap = 0; /* the number of Lyapunov numbers to be
computed <= vec_dim */
lyapzero = 2; /* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */

X_upper = pi; /* x scale */
X_lower = -pi;
Y_upper = pi; /* y scale */
Y_lower = -pi;
C1 = 1;
rho = 1.;
map = rotor; /* 31 */
}

LaserMap() { /* Ring Laser map:(Z= x+iy) Z -> rho+ C2*Z*exp{
i[C1 - C3/(1 + |Z|**2)] } */
double y_temp[6];
int base,
i;
double cs,

sn,
temp,
xn,
ynval, /* YN is a global variable */
C3over,
real,
imag;
double tempprime,
dfdx,
dfdy,
dgdx,
dgdy; /* first derivatives */
xn = y[0];
ynval = y[1];
C3over = C3 / (1 + xn * xn + ynval * ynval);
temp = C1 - C3over;
cs = C2 * cos(temp);
sn = C2 * sin(temp);
real = xn * cs - ynval * sn;
imag = xn * sn + ynval * cs;
y_temp[0] = rho + real;
y_temp[1] = imag;

if(num_lyap > 0) {
tempprime = 2 * C3over * C3over / C3;
dfdx = cs - imag * tempprime * xn;
dfdy = -sn - imag * tempprime * ynval;
dgdx = sn + real * tempprime * xn;
dgdy = cs + real * tempprime * ynval;

for(i = 0; i < num_lyap; i++) {
base = lyapzero + vec_dim * i;
y_temp[base]
= dfdx * y[base]
+ dfdy * y[base + 1];
y_temp[base + 1]
= dgdx * y[base]
+ dgdy * y[base + 1];
}
}

for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
y[i] = y_temp[i];
}

LaserInverseMap() { /* Ring Laser map:(Z= x+iy) Z -> rho+ C2*Z*exp{
i[C1 - C3/(1 + |Z|**2)] } */
double C2over,
psi,
ypsi,
xpsi,
dxpsi,
dypsi;
double arg,
dxarg,
dyarg;
double y_temp[6];
int base,
i;
double cs,
sn,
xn,
ynval;
double dfdx,
dfdy,
dgdx,
dgdy; /* first derivatives */


if(C2 == 0) {
PRINT "IMPOSSIBLE VALUE: C2 = 0; no inverse exists\n");
return;
}
xn = y[0];
ynval = y[1];
C2over = 1 / C2;
xpsi = (xn - rho) * C2over;
ypsi = ynval * C2over;
psi = 1 + xpsi * xpsi + ypsi * ypsi;
arg = C1 - C3 / psi;
cs = C2over * cos(arg);
sn = C2over * sin(arg);
y_temp[0] = (xn - rho) * cs - ynval * sn;
y_temp[1] = ynval * cs + (xn - rho) * sn;



if(num_lyap > 0) {
dxpsi = 2 * (xn - rho) * C2over;
dypsi = 2 * ynval * C2over;
dxarg = -C3 * dxpsi / (psi * psi);
dyarg = -C3 * dypsi / (psi * psi);
dfdx = cs - ((xn - rho) * sn + ynval * cs) * dxarg;
dfdy = ((-xn + rho) * sn - ynval * cs) * dyarg - cs;
dgdx = (-ynval * sn + (xn - rho) * cs) * dxarg - sn;
dgdy = cs + (-ynval * sn + (xn - rho) * cs) * dyarg;

for(i = 0; i < num_lyap; i++) {
base = lyapzero + vec_dim * i;
y_temp[base]
= dfdx * y[base]
+ dfdy * y[base + 1];
y_temp[base + 1]
= dgdx * y[base]
+ dgdy * y[base + 1];
}
}

for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
y[i] = y_temp[i];
}

initLaser() {
vec_dim = 2; /* the dimension of the Lyapunov vectors =
phase space dim */
num_lyap = 0; /* the number of Lyapunov numbers to be
computed <= vec_dim */
lyapzero = 2; /* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */

X_upper = 3; /* x scale */
X_lower = -1;
Y_upper = 4.5; /* y scale */
Y_lower = -4.0;
C2 = 0.9;
rho = 1.;
C3 = 6.;
C1 =.4;
map = LaserMap; /* 50 */
}

bob() { /* this routine defines the map */
int rand(); /* a Microsoft random integer generator */
double temp0,
temp1,
c,
s,
angle;

if(rand () < C2 * 32768.) {
/* rand() returns a random integer between 0
and 32768, so C2 should be a number between
0 and 1, so that the first process will be
carried out C2 of the time */
y[0] = 2.- y[0];
y[1] = -y[1];
}
else {
angle = C1 * twopi / 360.;
c = cos(angle);
s = sin(angle);
temp0 = rho * (c * y[0] + s * y[1]);
temp1 = rho * (-s * y[0] + c * y[1]);
y[0] = temp0;
y[1] = temp1;
}
}

init_bob() { /* this process initializes the constants and
with the first line tells the program the
name of the routine that defines the map;
since C1, C2, and rho are used in bob(),
htey must be initiallized in this routine;
if they are not, the program will not let be
accessed interactively and they will be
given default values, namely -9999. */
map = bob; /* this tells the program the name of the
routine that is evaluated; map is a pointer
to routines and bob() must be defined in the
file prior to this line; */
rho =.65;
C1 = 135.;
C2 =.5;
/* the following 4 determine the x and y scale */
X_lower = -3.;
X_upper = 5.; /* x scale */
Y_lower = -3.;
Y_upper = 3.; /* y scale */
}