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

/*************************** YBasins.c *************************************/
/********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
/********** ROUTINES in file:
SetNum_tests() sets num_tests which is the number of points to be checked by
checkTrajectory; we will not check points if a coordinate (the
zeroth coordinate) is -9999; we want num_tests-1 to be the
last storage vector for which that coordinate is not -9999
setframe(lower,upper,pedges) given two numbers, lower and upper, a new
interval (larger than the original if diameters[ScrnSec] > 0 )
is found and the ends are stored in pedges[0] and pedges[1]
where pedges[0] < pedges[1]
double lower,upper,*pedges;
BasinBoundary() initial points are plotted when storage vector 2 is not
set and the trajectory goes to infinity , that is when
(tau == 1 && y[eqn0 + eqn1 + zeroth] == -9999)
OR when the trajecory approaches storage vector 21
int checkTrajectory(EnterLevel,pVec)
This tests where the trajectory-that-starts-at-y is
initially and where it is after
each call of (*ChkTrajIteratee)(); checking the trajectory
for a maximum of "MaxChecks" times. The routine leaves 'time'
equal to the exit time (more or less). The routine:

returns 0 if it has gone thru the full set of MaxChecks checks
without getting near any of the specified attractors
and without leaving the box
returns 1 if it is ever outside the box frame[4];
returns k (2 <= k <= 7) if the trajectory of the point
approaches close to point yk[] /
int EnterLevel;
double *pVec;
int check_coord(pleft,Y)
pleft is a pointer to the value of the coordinate for
the left side of the box. We assume pleft+1 points to the
value of the right side. this function checks the two sides of
the box for a given number Y. This function returns
IN = 1 if point y is in interval,
OUT = 0 if outside, and
NOTSET = -9999 if neither box coord is defined. /
double *pleft,Y;
************************* STRADDLE ORBIT ROUTINES ************************
initSadStraddle()/* first we find out if ya and yb are set
and if so, what checkTrajectory returns for ya and yb;
the returned values are retA and retB;
if retA = retB, or if either equals 0, quit by decreasing
level; /

SetBasin() sets the vector basin[8] so that each coordinate is either
1 or 2 thereby indicating which attactor each yx[] point
is in; hence when getting a number that is returned from
checkTrajectory, the program can tell which basin it is in /
initBasStraddle() first we find out if ya and yb are set
and if some numbered storage vector >= 2 is set
and if so, what checkTrajectory returns for ya and yb;
the returned values are retA and retB;
if retA = retB, or if either equals 0, quit by decreasing
level; /
DoAB_Exit() checks where ya and yb are going and complains if they go
nowhere /
int AreAB_set() if ya and yb are set 1 is returned; otherwise 0 /
initChTraj()
BasinStraddle()
this routine takes the pair of vectors ya,yb and refines the pair
until they are "IsClose"; then it takes a single iterate of
the refined ya and yb so that the pair advances one iterate;
the method of refinement depends on BST_switch:
BST_switch = 1 means a random point between
ya and yb is chosen; otherwise the midpoint is chosen /
AccessBS() AccessBasinStraddle: this routine
takes the pair of vectors ya,yb and refines the pair
until they are "IsClose"; then it takes a single iterate of
the refined ya and yb so that the pair advances one iterate;
it is refined by taking lots of points between ya and yb and
as soon as any of them are found to be in basinB, move yb to
that point. If none of them are found to be in basinB, choose
a point near the midpoint that is known to be a basinA point
and move ya to it /
SaddleStraddle() we search for a trajectory that stays away from the basins;
this routine takes the pair of vectors ya,yb and refines the pair
until they are "IsClose"; then it takes a single iterate of
the refined ya and yb so that the pair advances one iterate /
int refineSadStrad()
this routine tries to find a refinement, i.e. a new pair
ya, yb, (one could be unchanged) closer together than the
starting pair; if it is successful, the pair is changed and
1 is returned;
it returns 0 if it finds an orbit that never leaves the domain;
it returns -1 if it cannot find a refinement; notice that all
the values of fraction are diadic and so are computed exactly
OneIterate(Y) iterates the point at Y[] using y /
double *Y;
set_y_fraction(fraction)
double fraction;
GiveUp()
triple()
******************************************************************************/

/* "never exits" case must be reworked to exit smoothly from program */
#include "yinclud.h"

#define OUT 0
#define IN 1
#define NOTSET -9999


double frame[4];
int num_tests;
static int isYkSet[9];
static int time;
static int retA,retB;
static int basin[8];
static int basinA,basinB;
static int intrptFlag = NO;
static double totalRefinementCalls,totalTime;
int modNum = 14;


SetNum_tests()/* sets num_tests which is the number of points to be checked by
checkTrajectory; we will not check points if a coordinate (the
zeroth coordinate) is -9999; we want num_tests-1 to be the
last storage vector for which that coordinate is not -9999 */
{
int k;
for (num_tests = NUM_Y;
y[eqn0 +(num_tests - 2)*eqn1 + zeroth] == -9999 && num_tests >= 2;
num_tests--)
;
isYkSet[8] = NO;
for (k = 2; k <= 7; k++)
if (y[eqn0 + eqn1*(k-1) + zeroth] != -9999. )
isYkSet[k] = YES;
else
isYkSet[k] = NO;
}


setframe(lower,upper,pedges)/* given two numbers, lower and upper, a new
interval (larger than the original if diameters[ScrnSec] > 0 )
is found and the ends are stored in pedges[0] and pedges[1]
where pedges[0] < pedges[1] */
double lower,upper,*pedges;
{
double temp;

pedges[0] = lower - diameters[ScrnSec]*(upper - lower);
pedges[1] = upper + diameters[ScrnSec]*(upper - lower);

if (lower > upper) /* switch so pedges[0] will be less than pedges[1]*/
{
temp = pedges[0];
pedges[0] = pedges[1];
pedges[1] = temp;
}
}
extern long initialdot;/* see Y.C */
BasinBoundary()/* initial points are plotted when storage vector 2 is not
set and the trajectory goes to infinity , that is when
(tau == 1 && y[eqn0 + eqn1 + zeroth] == -9999)
OR when the trajecory approaches storage vector 21 */
{
long jy;
int tau, level0=level, oldColor = color;
int cheqFreq();
double x0 ,x1, y0, y1, alittle,timeofday();

X_step = ( X_upper - X_lower)/xpixels;
Y_step = ( Y_upper - Y_lower)/ypixels;

level++;
if (level == PROCESS)
boot_crt();/* sets boot to 0, and if boot is 0, it clears
p[][] and the screen. It sets the scale */
SetBasin();/* sets basin[] */
initChTraj();

if (VertLine > 0 && VertLine < xpixels && refreshFlag == YES)
;
else
VertLine = 0;
refreshFlag = NO;

if (TDFreq > 0)
PRINT
"WARNING: picture is stored to disk every %ld vertical line(s) of dots.\n"
,TDFreq);

time0 = timeofday();
dot = dotAtTime0 = ((long) VertLine)*1000;
initialdot = ((long) VertLine)*1000;


initTime();
for (; (level == level0+1) && VertLine < xpixels; VertLine++)
{
if(timeFlag == YES || TDFreq > 0)
if(checkFreq() == YES)
{ /* means terminate; checkFreq also handles
occasional transmissions of picture to disk
*/
PRINT "process termination: OUT OF TIME \n");
level = level0 - 1;
continue;
}

alittle = .000001;
x0 = X_lower + X_step*(VertLine + alittle);
x1 = X_lower + X_step*(VertLine + 1 + alittle);
dot = VertLine*1000 ;/* 3 least significant digits = height */
for (jy = 0; jy < ypixels && level == level0 + 1; jy++)
{
dot++;
y0 = Y_upper - Y_step*(jy + alittle);
y1 = Y_upper - Y_step*(jy + 1 + alittle);
if (dim > 2) setequal(0,1,eqn0);
y[X_coord] = x0;
y[Y_coord] = y0;

tau = checkTrajectory(level,y);

if ( basin[tau] == 1)
{
color = (time % modNum)+1;
_setcolor(color);
block_plot(x0,x1,y0,y1);
/* in high resolution modes it is
possible that for some vertical lines of pixels
x0 = x1 within the resolution of the screen so
nothing will be plotted on the screen; it will
however be plotting in p[][] */
}
#ifndef MAINFRAME
if (ChkKeyStatus() != 0 || cycle > 0)
Interrupt();
#endif
}
}
if (VertLine != xpixels)
VertLine--;/* if we terminate in the middle of a picture, and
we begin again later, we want to begin at the
beginning of a line */
else
{
VertLine = 0;/* no harm to reset it if the entire picture
has been redrawn */
}
if ((onprint>=1) && (level == level0+1))
pprint(onprint);/* this sends picture to B: disk if disk = 1 */
#ifndef MAINFRAME
if (level == level0+1)
{
cycle = 1;
Interrupt();
}
#endif
level = level0;
setequal(0,1,eqn0);/* reinitializes y[] = y1 */
if (level < PROCESS)
{
scr_clr();/* clears screen but not pic[]; it is a function
in desmets pcio.a */
MainMenu();
}
color = oldColor;
}

whenAndWhere()/* tells where the trajectory through y1 goes and how long
it takes to get there; for interrupt W and command W */
{
double yTemp[EQUATIONS];
int ret,checkTrajectory();
int oldNumLyap = num_lyap;

num_lyap = 0;
store(yTemp,y,lyapzero);/* this routine operates on y[] */

SetBasin();/* sets basin[] */
initChTraj();

ret = checkTrajectory(level,y+eqn0);/* test y1 */
store(y,yTemp,lyapzero);/* this routine operates on y[] */
num_lyap = oldNumLyap;

if (ret == 0 && printer > 0)
PRINT
"Trajectory through small cross iterated MC=%d times, end point stored in ye\n"
,MaxChecks);
if (ret == 1 && printer > 0)
PRINT
"Trajectory through small cross diverges toward infinity in %d iterates. \n"
,time);

if (ret >= 2 && ret <= 7 && printer > 0)
PRINT
"Trajectory through small cross came within RA=%3.3lf of y%d in %d iterates\n"
,ra[ret],ret,time);
}

int checkTrajectory(EnterLevel,pVec)
/* This tests where the trajectory-that-starts-at-y is
initially and where it is after
each call of (*ChkTrajIteratee)(); checking the trajectory
for a maximum of "MaxChecks" times. The routine leaves 'time'
equal to the exit time (more or less). The routine:

returns 0 if it has gone thru the full set of MaxChecks checks
without getting near any of the specified attractors
and without leaving the box
returns 1 if it is ever outside the box frame[4];
returns k (2 <= k <= 7) if the trajectory of the point
approaches within RA of the point yk[] */
int EnterLevel;
double *pVec;
{
int k;
double distance(),*p1,*p2;


if (pVec != y)
store(y,pVec,lyapzero);/* this routine operates on y[] */

p1 = y+zeroth;
p2 = y + eqn0 + zeroth;
for (time = 0; (time < MaxChecks) && (level == EnterLevel); time++)
{

/* We now handle the case where the point is "diameters" outside
the screen region */
if( check_coord(frame, y[X_coord]) == OUT ||
check_coord(frame+2,y[Y_coord]) == OUT )
{
return (1);
}

for (k = 2; k < num_tests; k++)
if (isYkSet[k] == YES )/* has yk[] been set? */
{
if (distance(p1, p2 + (k-1)*eqn1, vec_dim)
< ra[k])
{
return(k);/* Y is close to kth test point */
}
}
(*ChkTrajIteratee)();
}
if (time == MaxChecks && MaxChecks >= 100 && level > 1)
{
store(ye,y,lyapzero);
PRINT
"A trajectory has been iterated MC times; final point is stored in ye[] \n");
PRINT
"Either set some y2...y7 to ye or use command MC to set MC to be about 20\n"
);

}
return (0);/* one way to get here is for Interrupt() to decrease
level via interrupt ; another way is for the trajectory to be
iterated the full MaxCheck iterates without meeting any of
the tests; */
}


int check_coord(pleft,Y)
/* pleft is a pointer to the value of the coordinate for
the left side of the box. We assume pleft+1 points to the
value of the right side. this function checks the two sides of
the box for a given number Y. This function returns
IN = 1 if point y is in interval,
OUT = 0 if outside, and
NOTSET = -9999 if neither box coord is defined. */
double *pleft,Y;
{
int ret;

if((*pleft == NOTSET) && (*(pleft+1) == NOTSET))
ret = NOTSET;
else
{
ret = ((*pleft != NOTSET)&&(*pleft > Y) ? OUT : IN );
if (ret == IN)
ret = ((*(pleft+1) != NOTSET)&&(*(pleft+1)) < Y ? OUT : IN);
}
return (ret);
}

/************************* STRADDLE ORBIT ROUTINES ************************/

initSadStraddle()/* first we find out if ya and yb are set
and if so, what checkTrajectory returns for ya and yb;
the returned values are retA and retB;
if retA = retB, or if either equals 0, quit by decreasing
level; */
{
totalRefinementCalls = totalTime = 0;

if ((AreAB_set()) == 0) /* 0 means>= one is not set;
if 0, then level is downgraded */
return;
initChTraj();

DoAB_Exit();/* If checkTrajectory returns 0, then level is downgraded*/
}

SetBasin()/* sets the vector basin[8] so that each coordinate is either
1 or 2 thereby indicating which attactor each yx[] point
is in; hence when getting a number that is returned from
checkTrajectory, the program can tell which basin it is in */
{
/* the number ret that checkTrajectory returns is:
0 if the trajectory never exits,
1 if it goes to infinity, i.e., leaves the box;
k if it enters the ball centered at yk[];

now we say basin is #1 if ret = 2,3, or 4 and
is #2 if ret = 5,6, or 7 or 0;
ret = 1 (infinity case: if it is ever outside the box
frame[4]) is basin #1 IF y2[0] is default = -9999.;
otherwise ret = 1 is basin #2 */


basin[2] = basin[3] = basin[4] = 1;
basin[0] = basin[5] = basin[6] = basin[7] = 2;
if (y[eqn0+eqn1] == -9999.)
basin[1] = 1;
else
basin[1] = 2;
/*
PRINT
"basin[0,7] = %d %d %d %d %d %d %d %d ]\n"
,basin[0],basin[1],basin[2],basin[3],basin[4],basin[5],basin[6],basin[7]);
*/
}

initBasStraddle()/* first we find out if ya and yb are set
and if some numbered storage vector >= 2 is set
and if so, what checkTrajectory returns for ya and yb;
the returned values are retA and retB;
if retA = retB, or if either equals 0, quit by decreasing
level; */
{
int initLevel = level;

SetBasin();
if (AreAB_set() == 0) /* 0 means>= one is not set */
return;
initChTraj();
DoAB_Exit();/* defines retA and retB */
basinA = basin[retA];
basinB = basin[retB];

initChTraj();
if (num_tests <= 2)/* set in initChTraj() */
{
PRINT
"\n @#$%^&*+# THERE IS A PROBLEM *&%$#*|@+ \n\n");
PRINT
"You must first tell the program where at least one attractor is. You must set"
);
PRINT
"\n at least one of the vectors 2 thru %d.\n",NUM_Y);
level = initLevel-1;
return;
}

if (basin[retA] == basin[retB])/* are they in the same basin? */
{
PRINT
"retA=5d; retB=%d\n",retA,retB);
PRINT
"basin[retA]=5d; basin[retB]=%d\n",basin[retA],basin[retB]);

PRINT
"Verrry baaad: stored points a & b are in the same basin %d\n"
,basin[retA]);
PRINT
"This occurs sometimes when the radius RA is too big, giving disks \n"
);
PRINT
" that actually lie in two basins\n");
level = initLevel-1;
}
store(y,yb,lyapzero);/* set y to equal yb so that the first point
plotted will not be a bad point; do not change
high coordinates of ya needed for exponents */
}


DoAB_Exit()/* checks where ya and yb are going and complains if they go
nowhere */
{
int initLevel = level;

storeNumLyap = num_lyap;
num_lyap = 0;/* do not compute exponents during this phase */

retA = checkTrajectory(level,ya);
if (printer >= 2)
PRINT
"ya returns attractor# %d in time %d \n"
,retA,time);
retB = checkTrajectory(level,yb);
if (printer >= 2)
PRINT
"yb returns attractor# %d in time %d \n\n"
,retB,time);
num_lyap = storeNumLyap;

if (retA == 0)
{PRINT
"stored point a returns zero: its trajectory goes to no attractor \n");
level = initLevel - 1;
return;
}
if (retB == 0)
{PRINT
"stored point b returns zero: its trajectory goes to no attractor \n");
level = initLevel-1;
return;
}
}


int AreAB_set()/* if ya and yb are set 1 is returned; otherwise 0 */
{
if (ya[zeroth] == -9999)
{ PRINT
"stored point a is still at the default value; cannot proceed \n");
level--;
return(0);
}
if (yb[zeroth] == -9999)
{ PRINT
"stored point b is still at the default value; cannot proceed \n");
level--;
return(0);
}
return (1);
}

initChTraj()
{
int k;
/* checkTrajectory checks to see if the trajectory has left the box with
edges in frame[] so we have so set frame[] */
setframe(X_Lo[ScrnSec],X_Up[ScrnSec],frame);
setframe(Y_Lo[ScrnSec],Y_Up[ScrnSec],frame+2);/* the box with
coordinates
frame[0] < x < frame[1]; frame[2] < y < frame[3]
is larger than the screen box by a factor 1 + 2* diameters */
SetNum_tests();/* num_tests is the number of points to be checked by
checkTrajectory; we will not check points if a coordinate is -9999 */

for (k=2; k<9; k++)
ra[k] = rad_attr;
if (ra2 != -9999.)/* ra2 and ra5 are set using commands RA2 & RA5
in YCOMB.C; see also YPICKMAP.C*/
ra[2] = ra2;
if (ra5 != -9999.)
ra[5] = ra5;
}


BasinStraddle()
/* this routine takes the pair of vectors ya,yb and refines the pair
until they are "IsClose"; then it takes a single iterate of
the refined ya and yb so that the pair advances one iterate;
the method of refinement depends on BST_switch:
BST_switch = 1 means a random point between
ya and yb is chosen; otherwise the midpoint is chosen */
{
int random997(); /* in YTOOLS.C Returns an integer between 0 and 997 */
int ret;
int initLevel = level;
double dist,distance(),fraction;

storeNumLyap = num_lyap;
num_lyap = 0;

ret = retA; /* if dist below turns out to be small, ret needs a
non zero value */


/* now the procedure begins */
while ((dist = distance(ya+zeroth,yb+zeroth,vec_dim)) > IsClose
&& level == initLevel )
{
/* try to refine the pair because they are NOT IsClose */
/* SET y EQUAL TO THE AVERAGE OF ya AND yb */
if (BST_switch == 1)
fraction = (1.0 + random997())/999.;
else
fraction = 1./2.;

set_y_fraction(fraction);
if (ChkKeyStatus() != 0 )
{
Interrupt();
if (level < initLevel)
continue;/* process terminated via keyboard
Interrupt */
}

/* Now see where y goes */
ret = checkTrajectory(level,y);
scr_rowcol(0,0);
if (ret == 0)
{
erase_line();
PRINT
"\nret = 0: y didn't go to any attractor; program doesn't know what to do;\n");
PRINT
"return to single step mode... but we are still in BasinStraddle() \n");
cycle = 1; /* single step mode; */
Interrupt();
}
if (printer == 3)
{
PRINT
"\nab dist=%5.5le, time=%d, dot# =%ld ret=%d (to suppress printing, hit 2)\n",
dist,time,dot,ret);

}
/* now consider case where we fail to refine */
if (dist < 10.*IsClose && time == 0)
{
PRINT
"ya[] and yb[] are now in region %d so one has switched basins\n",ret);
GiveUp();
}



set_y_fraction(fraction);/* RECALL y */

/* REFINE ya yb line */
/* CHOOSE A NEW ya OR yb TO BE THIS MIDPOINT y */
if (basin[ret] == basinA)
store(ya,y,lyapzero);/* put y in ya; do not change
high coordinates of ya needed for exponents */
else
if (ret != 0 && basin[ret] == basinB)/* ret = 0 means the
trajectory did not go to a basin, possibly
because the number of iterates checked = MC
was too low; the routine should stop without
changing yb */
store(yb,y,lyapzero);

}
/* if (ret != 0) changed 5/14/87; this is why we gave a default value
to ret at top of this routine */

OneIterate(ya);
num_lyap = storeNumLyap;
OneIterate(yb);
/* Notice that yb[] = y[] now and this is the point that gets plotted and if
Lyapunov exponents are being calculated, they are essentially
the exponents of yb */
}

AccessBS()/* AccessBasinStraddle: this routine
takes the pair of vectors ya,yb and refines the pair
until they are "IsClose"; then it takes a single iterate of
the refined ya and yb so that the pair advances one iterate;
it is refined by taking lots of points between ya and yb and
as soon as any of them are found to be in basinB, move yb to
that point. If none of them are found to be in basinB, choose
a point near the midpoint that is known to be a basinA point
and move ya to it */
{
int initLevel = level;
int ret, N;
double dist,distance(),frac,GoldenMean;
GoldenMean = 2./(1+ sqrt(5.));/* note this is less than 1 */

storeNumLyap = num_lyap;
num_lyap = 0;
ret = retA; /* if dist below turns out to be small, ret needs a
non zero value */

/* now the procedure begins */
while ((dist = distance(ya+zeroth,yb+zeroth,vec_dim)) > IsClose
&& level == initLevel )
{
/* try to refine the pair because they are NOT IsClose */
/* SET y EQUAL TO THE AVERAGE OF ya AND yb */

for (frac = GoldenMean, N = 0; N < divisions && level == initLevel;
N++,frac = frac+GoldenMean)
{
while (frac >= 1)
frac = frac - 1;
set_y_fraction(frac);
/* Now see where y goes */
ret = checkTrajectory(level,y);
if (ret == 0)
{
scr_rowcol(0,0);
erase_line();
PRINT
"\nret = 0: y didn't go to any attractor; program doesn't know what to do;\n");
PRINT
"return to single step mode... but we are still in BasinStraddle() \n");
cycle = 1; /* single step mode; */
Interrupt();
}/* end if */
if (printer == 3)
{
scr_rowcol(0,0);
PRINT
"\nab dist=%5.5le, time=%d, dot# =%ld ret=%d (to suppress printing, hit 2)\n",
dist,time,dot,ret);
}/* end if */
/* now consider case where the intermediate point is close to
an attractor (i.e. time == 0 */
if (dist < 10.*IsClose && time == 0)
{
PRINT
"ya[] and yb[] are now in region %d so one has switched basins\n",ret);
GiveUp();
}/* end if */

/* REFINE ya yb line if basin[ret] == basinB */
if (basin[ret] == basinB)/* note this includes the case ret = 0
so it could cause problems */
{
set_y_fraction(frac);/* RECALL y */
store(yb,y,lyapzero);
break;/* terminates for() */
}/* end if */
}/* end for */
if (basin[ret] == basinA)
{
set_y_fraction(GoldenMean);/* RECALL a y that gave a
ret indicating basinA */
store(ya,y,lyapzero);
}
}/* end while */
/* if (ret != 0) changed 5/14/87; this is why we gave a default value
to ret at top of this routine */

OneIterate(ya);
num_lyap = storeNumLyap;
OneIterate(yb);
/* Notice that yb[] = y[] now and this is the point that gets plotted and if
Lyapunov exponents are being calculated, they are essentially
the exponents of yb */
}


SaddleStraddle()/* we search for a trajectory that stays away from the basins;
this routine takes the pair of vectors ya,yb and refines the pair
until they are "IsClose"; then it takes a single iterate of
the refined ya and yb so that the pair advances one iterate */
{
int initLevel = level;
int success;
char *comment;
double distance(),RefCallsPerDot,TimePerPlot,dist;

storeNumLyap = num_lyap;
num_lyap = 0;
success = 1;/* if ya and yb are close, skip stuff in while and
iterate again */
while ((dist = distance(ya+zeroth,yb+zeroth,vec_dim)) > IsClose
&& level == initLevel )
{
/*try to refine the pair because they are NOT within IsClose */
success = refineSadStrad();

if (success == 1)
comment = "success";
if (success == 2)
comment = "KEYBOARD TERMINATION";
if (success == 0)
comment = "Never exits";
if (success == -1)
comment = "refinement failure";
if (printer >= 2 )
{
dist = distance(ya+zeroth,yb+zeroth,vec_dim);
scr_rowcol(0,0);
erase_line();
PRINT
"ab dist=%5.5le, time=%d, dot# =%ld, %s (to suppress printing, hit 2)\n",
dist,MaxTime,dot,comment);
if (dot > 0)
{
RefCallsPerDot = totalRefinementCalls/dot;
TimePerPlot = totalTime/dot;
erase_line();
PRINT
"refinement calls per dot = %lf;av. time per plot=%lf;\n"
,RefCallsPerDot,TimePerPlot);
erase_line();
PRINT
"totalRefinementCalls=%3.3lf;totalTime=%3.3lf; \n"
,totalRefinementCalls,totalTime);
}
}
/* now consider case where we fail to refine */
if (dist < 10.*IsClose && MaxTime == 0)
{
PRINT
"ya[] and yb[] have MaxTime = 0 so process should not proceed\n");
PRINT
"This suggests points in one ball can actually go to another region\n");
GiveUp();/* return to single stepping */
}
}

if (success == 1)/* this is the good case */
{
totalTime = totalTime + time;
OneIterate(ya);
num_lyap = storeNumLyap;
OneIterate(yb);
}
num_lyap = storeNumLyap;
/* Notice that yb[] = y[] now and this is the point that gets plotted and if
Lyapunov exponents are being calculated, they are essentially
the exponents of yb */
}



int refineSadStrad()
/* this routine tries to find a refinement, i.e. a new pair
ya, yb, (one could be unchanged) closer together than the
starting pair; if it is successful, the pair is changed and
1 is returned;
it returns 0 if it finds an orbit that never leaves the domain;
it returns -1 if it cannot find a refinement; notice that all
the values of fraction are diadic and so are computed exactly*/
{
double oldNumLyap;
double slope, slope0;/* for access. basin strad. */
double Ytemp[EQUATIONS],fraction=0;
int initLevel = level;
int FindTime();
int numA,numB,monotone,decMonotone,testA;
int ndivisions; /* a working multiple of divisions */
int times[1001],N,firstMax;/* firstMax records the first time out of
NumSegment+1 times that the value MaxTime is taken on */

totalRefinementCalls++;/* counts the number of refinements */
for (ndivisions = divisions; ndivisions < 1000
&& level == initLevel ;
ndivisions *= 2)/* integer 1000 can't exceed dim times */
{
scr_rowcol(6,0);
MaxTime = 0;
decMonotone = monotone = 0;
firstMax = 0;
for (N = 0; N <= ndivisions && level == initLevel ; N++)
{
#ifndef MAINFRAME
if (ChkKeyStatus() != 0 )
#endif /* ifndef MAINFRAME */
{
set_y_fraction (fraction);
oldNumLyap = num_lyap;
Interrupt();/* this may reset the number of Lyapunov
exponens via a #L command; thus we must
check to see if num_lyap has been changed;
if so, storeNumLyap is what should have
been set, since num_lyap is supposed to be
0 at this point in the program */
if (oldNumLyap != num_lyap)
{
storeNumLyap = num_lyap;
num_lyap = 0;
}
if (level < initLevel)
return (2) ; /* process terminated via keyboard
interrupt */
}
fraction = ((double)N)/ndivisions;/* we want a double
precision number */
FindTime(fraction);/* ret is what checkTra.. returns;
"time" is set by FindTime() */
times[N] = time;
if (decMonotone == N-1 && times[N-1] >= times[N])
{
decMonotone++;/* measures how long the sequence is
non-increasing */
monotone = decMonotone;/* monotone is supposed to
measure how long times[] is non-decreasing
AFTER the initial period of decrease */
}
if (monotone == N-1 && times[N-1] <= times[N])
monotone++;/* monotone is supposed to
measure how long times[] is non-decreasing
AFTER the initial period of decrease */


if (printer >= 3 || gameFlag)
{
PRINT "(%d,%d) ",N,time);
}
if (time > MaxTime)
{
firstMax = N;
MaxTime = time;
}
}/* end inner for */
if (printer >= 2 || gameFlag == YES)
{
PRINT
"\n max point = (%d,%d)\n"
,firstMax,MaxTime);
if (decMonotone == N-1 && times[0] > times[decMonotone])
PRINT
"monotone decreasing thru (%d,%d); then \n",
decMonotone,times[decMonotone]);
else
PRINT
"constant thru (%d,%d); then \n",
decMonotone,times[decMonotone]);
PRINT
"monotone increasing thru (%d,%d)\n"
,monotone,times[monotone]);
}
#ifndef MAINFRAME /* cant play interactive game in batch mode */
if (gameFlag == YES)
{
numA = numB = 0;

while (numA == numB)
{
PRINT
"Input two integers >= 0 and <= %d for new ya and yb; NO COMMA between them \n"
,divisions);
testForInterrupt();
if (level == initLevel)
{
if (intrptFlag == YES)
PRINT
"Input two integers >= 0 and <= %d for new ya and yb; NO COMMA between them \n"
,divisions);
scanf("%d %d",&numA,&numB);
}
else
break;
boot = 1;

boot_crt();
}
}/* end gameFlag if */
#endif /* ifndef MAINFRAME */

/* ASST: */
else if (accessFlag == YES)/* decMonotone reports the last in a
series of times[] values that do not strictly
increase; monotone
reports the last subscript in the following series
of non decreasing times[] values; hence
monotone > decMonotone;
*/
{
numA = decMonotone;
numB = monotone+1;
if (monotone > ndivisions)
numB = ndivisions;
testA = ndivisions/3;
if (testA < monotone && testA > numA)
{
numA = testA;
if (printer == 3)
PRINT
"choice of numA = ndiv/3 = %d; decMon=%d;\n",numA,decMonotone);
}
if (numA == 0 && numB == ndivisions)/* i.e., no change, then
we try a refinement, but it looks bad */
break;/* break out of for */
PRINT "Interval refined from (0,%d) to (%d,%d); dot = %ld;\n"
,ndivisions,numA,numB,dot);
/*******/
}
else /* SADDLE STRADDLE TRAJECTORY;
in this case, we find the point in the grid where the time
is maximized; then we find the grid points, one on each side
which maximize the slope of the times[] function; these grid
points are numA and numB, which lets us define ya and yb */
{
numA = 0;
slope0 = 0;
if (firstMax > 0)
for (N = 0; N < firstMax; N++)
{
slope = (MaxTime-times[N])/(firstMax-N);
if (slope >= slope0)
{
slope0 = slope;
numA = N;
}
}
else
numA = -divisions;/* this causes the interval to
expand */
numB = divisions;
slope0 = 0;
if (firstMax < divisions)
{
for (N = divisions; N > firstMax; N--)
{
slope = (MaxTime-times[N])/(N - firstMax);
if (slope >= slope0)
{
slope0 = slope;
numB = N;
}
}
if (slope == 0)
numB = 2*divisions;/* this causes the interval
to expand */
}
else
numB = 2*divisions;/* this causes the interval to
expand */
}
if (level == initLevel)
if (printer >= 2)
PRINT "Selected grid points are: %d and %d \n",numA,numB);

/* numA and numB are now set so ya and yb are now changed */
fraction = ((double)numB)/ndivisions; /*we want a double
precision number */
set_y_fraction (fraction);
store(Ytemp,y,lyapzero);
fraction = ((double)numA)/ndivisions;
set_y_fraction (fraction);
store(ya,y,lyapzero);
store(yb,Ytemp,lyapzero);
return (1);/* means refinement is achieved; */
} /* end outer for */
return (1);/* means refinement is achieved; */
}






OneIterate(Y)/* iterates the point at Y[] using y */
double *Y;
{
store(y,Y,eqn0);
(*StraddleIteratee)();
store(Y,y,eqn0);
}

set_y_fraction(fraction)
double fraction;
{
int k;
for (k = 0; k < lyapzero; k++)
y[k] = (1-fraction) * ya[k] + fraction * yb[k];
}

GiveUp()
{
erase_line();
PRINT
"Program now enters SINGLE STEP because of this failure.\n");
if (storeNumLyap != 0)
num_lyap = storeNumLyap;

cycle = 1; /* single step mode */
Interrupt();
}

int FindTime(fraction)
double fraction;
{
int ret;

set_y_fraction(fraction);
ret = checkTrajectory(level,y);
if (ret == 0) /* == 0 means trajectory did not exit */
{
set_y_fraction(fraction);
/* this will be the value of y that is in y
when SaddleStraddle() dies */
}

return (ret);/* == 0 if point does not go anywhere */
}

triple()
{
double Y[EQUATIONS];

set_y_fraction(-1.);
store(Y,y,eqn0);
set_y_fraction(2.);
store(yb,y,lyapzero);
store(ya,Y,lyapzero);
}

#ifndef MAINFRAME
testForInterrupt()
/* For detecting non digit when a digit is expected;
any character except \n is put back in the stream "input";
for anything other than a digit, it returns YES meaning YES
otherwise NO;*/
{
int ch;
int ChkKeyStatus();
while (ChkKeyStatus() == 0)
;
while ((ch = ci())== '\n')/* there may be a return still in
the pipeline from entering the command so this one should
not cause an abort */
;
ChkKeyStatus();
if (ch >= '0' && ch <= '9')
{
intrptFlag = NO;
ungetc(ch,stdin);
return;
}
else
{
intrptFlag = YES;
PRINT "PAUSE MODE; interrupt is now expected\n");
cycle = 1;
Interrupt();
return;
}
}
#endif /* ifndef MAINFRAME */