Category : C Source Code
Archive   : DYNAM_S2.ZIP
Filename : YU.C
/********************************** YU.C ***********************************/
#include "yinclud.h"
static double diagonal; /* for unstable manifold routines in this file
*/
/* int goodIterate; the number of iterates before leaving the box
(set in terms of diameters[] ) */
static double ddx,
ddy;
static double PreviousAdd,
oneOverX,
oneOverY;
static int multiplier;
static double jump,
oldX,
oldY,
dif;
static int itermax;
LineIterator(iterMax, iterMultiplier)/* The images of the line from
x1_line,y1_line to x2_line,y2_line are
drawn, where iterMultiplier is a number >=1
such that the line is iterated in multiples
of iterMultiplier; in the computation of
unstable manifolds of periodic orbits, it
should be 2*period where period is the
period of the orbit; the factor of two is
because the orbit might have a negative
expanding eigenvalue; iterMax is the maximum
number of times "iter" is increased; the
actual maximum number of iterates is
iterMax*iterMultiplier; routines in this
file: set_line(),LineForB(), and LineForU()
set the points(line_xi,line_yi) for i = 1
and 2. */
int iterMax,
iterMultiplier;
{
int isInBox(), level0;
int checkFreq();
double minAdd = 1;
double getAdd();
double fracY[EQUATIONS];
double oldLength = 0.,
entropy,
difX,
difY;
/* the following are for restarting this routine */
/* dot is just for keeping track of how much has been done */
if(refreshFlag == YES) {
if(dot > 0 && dot < dots)
/* dot >0 is most likely to occur when using
command "cont" */
;
}
else {
dot = 0;
iter = 1;
}
/* dot is just for keeping track of how much has been done */
initTime();
itermax = iterMax;
multiplier = iterMultiplier;/* we need this number globally */
level0 = ++level;
if(level == PROCESS && refreshFlag == NO)
boot_crt();
store(y, y + eqn0, eqn1);
oldX = y[X_coord];
oldY = y[Y_coord];
jump = USTEP; /* this is a fraction of the screen's diagonal
length(currently .001) indicating how big a
jump is allowed so that we do not move more
than 1 pixel */
difX = line_x1 - line_x2;
difY = line_y1 - line_y2;
length = sqrt(difX * difX + difY * difY);
/* this will tally the total length plotted so
far; initially it is the length of the
initial line segment */
diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
+ (Y_upper - Y_lower) * (Y_upper - Y_lower));
boxConstants(); /* this defines constants needed by isInBox()
for checking if y[] is in the box */
oldGoodIt = 0;
goodIterate = 1;
iter = 1;
add = USTEP;
store(fracY, y, dim); /* needed for getAdd() */
add = getAdd(add, fracY);/* one iterate; only place called */
for(; level == level0 && iter <= itermax; iter++) {
goodIterate = iter;
/* first we compute the initial jump size "add" */
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;
}
scr_rowcol(4, 0);
if(printer >= 2) {
if(iter > 1 && length > 0.&& oldLength > 0.) {
entropy = log(length / oldLength);
"length(%d)=length of %dth image of segment=%lf \n",
iter - 1, iter - 1, length * diagonal);
"Crude Estimate of topological entropy = log(length(%d)/length(%d))=%5.5lf \n"
,iter - 1, iter - 2, entropy);
}
else
PRINT "length(0)=%lf \n", length);
}
oldLength = length;
length = 0.;
/* MAIN INNER LOOP */
if(refreshFlag == YES)
refreshFlag = NO;
else
frac = 0;/* this occurs always except possibly once */
for(; level == level0 && frac < 1.0;
frac = ((frac + add >= 1.0) ? 1.0 : frac + add)) {
checkForInterrupt();
dot++; /* can be deleted; just for keeping track of
how much has been done and how fast it is */
/* when frac == 0, goodIterate and oldGoodIterate will not be
equal since goodIterate has just been increased beyond
previous values */
getGoodIt();
/* this is the only place in the "for(frac..."
loop that goodIterate is set; this also sets
y according to frac and iterates as many
times up to iter as possible while staying
in the box; at the end, y[] is this final
image point that is in the box and
goodIterate is the number of good iterates
*/
while(goodIterate >= oldGoodIt + 2) {
scr_rowcol(7, 0);
"good Iterate just jumped by 2\n");
frac = frac - add;
add = add / 5.;
frac = frac + add;
getGoodIt();
}
store(fracY, y, dim);/* fracY is used by revise... */
if (printer >= 3)
{
scr_rowcol(9,0);
"%d good; iter=%d frac=%lf,add=%le \n",goodIterate,iter,frac,add);
}
if (printer >= 1 && minAdd > add)
{
scr_rowcol(8,0);
"min add so far =%le \n",add);
}
minAdd = (minAdd > add ? add : minAdd ) ;
reviseAddAndJump(goodIterate, fracY);
/* if(goodIterate <= iter-1) we test to see if frac+add gives
a point that can be iterated one more time; if it can be,
then we will be decrease add by implementing the instruction
: reviseAddAndJump(goodIterate+1,fracY) which should produce
a smaller "add" and just to be safer, we "add = add/3.;" */
if(goodIterate <= iter - 1) {
process(1);
/* one more iterate from y[] as finally set by
revise...() */
if(isInBox (y[X_coord], y[Y_coord]) == YES) {
store(y, fracY, dim);
process(1);
/* one more iterate from y[] as finally set by
GetGoodIt() */
store(fracY, y, dim);
/* now fracY is as set by GetGoodIt() PLUS 1
iterate */
reviseAddAndJump(goodIterate + 1, fracY);
add = add / 3.;
}
}
if(goodIterate == iter || oldGoodIt == iter) {
store(y, fracY, dim);
plotSeg();
length = length + dif;
}
oldGoodIt = goodIterate;
/* the following two coordinates result from iterating y[] iter
times PROVIDED goodIterate = either iter or iter -1; then
they can be used for plotting */
oldX = y[X_coord];
oldY = y[Y_coord];
}
}
level = level0 - 1;
if(level < PROCESS)
MainMenu();
}
double jumpSize() { /* returns the maximum of the horizontal and
vertical fractions of the screen between y
and(oldX,oldY) */
double fabs(), difX, difY;
difX = fabs((y[X_coord] - oldX) * oneOverX);
difY = fabs((y[Y_coord] - oldY) * oneOverY);
if(difX > difY)
return(difX);
return(difY);
}
reviseAddAndJump(its, pY) /* given an estimate "add" of how much frac
should be changed by, try it out using
MakeAJump() which tells us what fraction of
the screen we actually move; if that move is
too big, (that is, if the updated estimate
is less than .8 of add), then do it again
using the new add; the desired jump is
jump*numPixels so we change add by the
factor jump*numPixels/dif. ALERT: if
goodIterate < its at any point, we stop
cycling */
int its;
double *pY;
{
double MakeAJump();
PreviousAdd = add;
for(;;) { /* change "add"; if it decreases, don't
decrease it too much in one step; but keep
re-evaluating while it decreases; decreasing
too much causes a problem when two points
differ primarily because they are being
calculated mod something */
dif = MakeAJump(its, frac + add, pY);
add =.2 * add +.8 * add * numPixels * jump / dif;
/* if add is decreased by more than 20%, go through another time to
refine add */
if(add <= 0.5 * PreviousAdd) {
add = 0.5 * PreviousAdd;
PreviousAdd = add;
}
else
break;
/* checkForInterrupt(); seems possibly unnecessary */
}
if(add > 2.* PreviousAdd) {/* this multiplier should be at least the
reciprocal of multiplier of PreviousAdd */
add = 2.* PreviousAdd;/* don't increase add too fast */
}
}
double MakeAJump(times, fraction, pY)
/* The process is iterated "times" times and
compares how far the new y[] is in the plane
from the vector pY[], and returns that
distance(unless the distance is too small;
then it returns .00...01 ); */
int times;
double fraction,
*pY;
{
double difX,
difY,
ret;
setY(fraction);
process(times);
difX = y[X_coord] - pY[X_coord];
difY = y[Y_coord] - pY[Y_coord];
ret = sqrt(difX * difX + difY * difY) / diagonal;
if(ret <.000001)
ret =.000001;
return(ret);
}
process(times) /* the inner loop (which itself iterates
multipier times) is iterated up to "iter"
times and in fact = "iter" times unless the
trajectory leaves the box, in which case
iterating is terminated. The position y[] at
exit time is its value after "times"
iterates */
int times;
{
int i,
j;
for(i = 1; i <= times; i++) {
for(j = 0; j < multiplier; j++)
(*iteratee) ();
}
}
int isInBox(xval, yval) /* This routine uses the same arithmetic
procedures as plot to see if
(y[X_coord],y[Y_coord]) is within the window
of ScreenSection and if "diameter" is >= 0
and it is not within the window, the routine
checks to see if it is within
"diameters[ScrnSec]" times the width and
within "diameters[ScrnSec] times the height
of the screen. If it is not close enough to
the screen based on this criterion, the
program returns NO, otherwise YES */
double xval,
yval;
{
if((X_upper - xval) * (X_lower - xval) > ddx
|| (Y_upper - yval) * (Y_lower - yval) > ddy) {
return(NO);
}
return(YES);
}
int isInCRT(xval, yval) /* This routine uses the same arithmetic
procedures as plot to see if
(y[X_coord],y[Y_coord]) is within the window
of ScreenSection. If it is not in the
screen, the program returns NO, otherwise
YES */
double xval,
yval;
{
if((X_upper - xval) * (X_lower - xval) > 0
|| (Y_upper - yval) * (Y_lower - yval) > 0) {
return(NO);
}
return(YES);
}
boxConstants() { /* this defines constants needed by isInBox()
for checking if((X_upper-xval)*(X_lower-xval) >
ddx || (Y_upper-yval)*(Y_lower-yval) > ddy )
the constants ddx and ddy must be reset if
the scale or box is changed */
double d,
dd,
xdif,
ydif;
xdif = X_upper - X_lower;
ydif = Y_upper - Y_lower;
oneOverX = 1 / (X_upper - X_lower);
oneOverY = 1 / (Y_upper - Y_lower);
d = diameters[ScrnSec];
dd = d * (1 + d);
ddx = dd * xdif * xdif;
ddy = dd * ydif * ydif;
}
checkForInterrupt() {
#ifndef MAINFRAME
if(ChkKeyStatus () != 0 || cycle > 0) {
Interrupt();
boxConstants();/* this defines constants needed by isInBox()
for checking if y[] is in the box, and
Interrupt can change the required constants
*/
}
#endif
}
plotSeg() { /* plots one or more points depending on
numPixels and modFlag */
double Xp0,
Yp0;
double jumpSize(), js;
int isInCRT(), nPix = numPixels;
int i;
if(goodIterate == iter - 1 || oldGoodIt == iter - 1)
nPix = 10 * nPix;
js = jumpSize();
if(js >.05 + numPixels *.001)
return;
if(cross0flag == 1) { /* this is for preparing to plot a cross */
plotX = y[X_coord];
plotY = y[Y_coord];
}
if(nPix > 1) { /* this segment can screw up if either or both
points were altered by moduloAB() */
Xp0 = y[X_coord];
Yp0 = y[Y_coord];
/* modFlag == YES means that a number has been changed by the
moduloAB() */
for(i = 1; i <= nPix; i++) {
modFlag = NO;
plotX = (Xp0 * i + oldX * (nPix - i)) / nPix;
plotY = (Yp0 * i + oldY * (nPix - i)) / nPix;
if(isInCRT (plotX, plotY) == YES)
plot(plotX, plotY);
}
}
else { /* if numPixels is not greater than 1 */
modFlag = NO;
if(isInCRT (y[X_coord], y[Y_coord]) == YES)
plot(y[X_coord], y[Y_coord]);
}
}
getGoodIt() { /* sets y according to frac and iterates as
many times up to iter as possible while
staying in the box; at the end, y[] is this
final image point and goodIterate is the
number of iterates */
double goodY[EQUATIONS];
int isInBox();
int i,
j;
goodIterate = 0; /* this was = 1 in a buggy earlier version
9/7/88 */
setY(frac);
store(goodY, y, dim);
for(i = 1; i <= iter; i++) {
for(j = 0; j < multiplier; j++)
(*iteratee) ();
if(isInBox (y[X_coord], y[Y_coord]) == YES) {
store(goodY, y, dim);
goodIterate = i;
}
else {
store(y, goodY, dim);
return;
}
}
}
setY(fraction)
double fraction;
{
int i;
for(i = zeroth; i < zeroth + vec_dim; i++)
y[i] = ya[i] + fraction * (yb[i] - ya[i]);
/*** y[X_coord] = line_x1 + fraction*(line_x2 - line_x1);
y[Y_coord] = line_y1 + fraction*(line_y2 - line_y1); ***/
}
double getAdd(add, pY) /* returns estimate for add */
double add,
*pY;
{
int n;
double MakeAJump();
getGoodIt();
/* now estimate add and then refine the estimate */
for(n = 0; n < 2; n++) {
dif = MakeAJump(goodIterate, frac + add, pY);
/* should be about equal to jump*numPixels */
add = add * numPixels * USTEP * diagonal / dif;
}
return(add);
}
#ifdef OBSOLETE
set_line(xshift1, xshift2, yshift1, yshift2)
/* used for drawing lines and unstable
manifolds */
double xshift1,
xshift2,
yshift1,
yshift2;
{
line_x1 = y[eqn0 + X_coord] + xshift1;
line_x2 = y[eqn0 + X_coord] + xshift2;
line_y1 = y[eqn0 + Y_coord] + yshift1;
line_y2 = y[eqn0 + Y_coord] + yshift2;
}
#endif
LineForU(iterMultiplier) /* this routine sets the ends of a line from ya
to yb for computing an unstable manifold
assuming that period NP in Newton's Method
divides the period of the orbit and y1[] is
at a point of the orbit */
int iterMultiplier;
{
double diagonal,
distance(), ratio;
double yy[EQUATIONS],
difX,
difY,
dif;
int i,
J;
diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
+ (Y_upper - Y_lower) * (Y_upper - Y_lower));
/*******
if(printer >= 3)
"initial:y1(x) = %11.11lf, y() = %11.11lf, dif=%11.11lf\n"
,y[eqn0+X_coord],y[X_coord],y[X_coord]-y[eqn0+X_coord]);
*********/
/* the following is to set the time coordinates of ya and yb so they will be
equal; other coordinates copied are now irrelevant */
store(ya, y, zeroth);
store(yb, y, zeroth);
/* iterate until y is more than dif/(*USTEP*diagonal) from y1 */
J = 0; /* counts iterates to prevent loop from never
ending */
do {
store(yy, y, dim);
for(i = 0; i < iterMultiplier; i++)
(*iteratee) ();
difX = y[X_coord] - y[eqn0 + X_coord];
difY = y[Y_coord] - y[eqn0 + Y_coord];
dif = sqrt(difX * difX + difY * difY);
ratio = dif / (numPixels * USTEP * diagonal);
/* desired ratio = 1 */
if(printer >= 1 && J == 0 && ratio > 5.)
/* ratio should start out small */
warning1 ();
/********* if(printer >= 3)
{
"J = %d: y1(x) = %11.11lf, y() = %11.11lf, dif=%11.11lf\n"
,J,y[eqn0+X_coord],y[X_coord],
y[X_coord]-y[eqn0+X_coord]);
"iteration %d: length = %lf length/desired len.= %lf difX=%lf \n"
,J,dif,ratio,difX);
} *********/
} while(++J < 100 && ratio < 1.);
if(J == 100)
"Something is wrong. The point may not be unstable \n");
if(printer >= 2)
PRINT "diagonal = %lf ratio = %lf \n", diagonal, ratio);
/* now rescale yy so that is iterates to a point dif/(USTEP*diagonal) from y1*/
if(printer >= 3)
PRINT "before first rescaling, ratio is now %lf\n", ratio);
if(ratio > 0) {
for(i = zeroth; i < lyapzero; i++)
yy[i] = y[eqn0 + i] + (yy[i] - y[eqn0 + i]) / ratio;
store(y, yy, dim);
for(i = 0; i < iterMultiplier; i++)
(*iteratee) ();
difX = y[X_coord] - y[eqn0 + X_coord];
difY = y[Y_coord] - y[eqn0 + Y_coord];
}
dif = sqrt(difX * difX + difY * difY);
ratio = dif / (numPixels * USTEP * diagonal);/* desired ratio = 1 */
if(printer >= 3)
PRINT "after first rescaling, ratio is now %lf\n", ratio);
/* now rescale again in an effort to get closer to the desired length */
for(i = zeroth; i < lyapzero; i++)
yy[i] = y[eqn0 + i] + (yy[i] - y[eqn0 + i]) / ratio;
store(y, yy, dim);
for(i = 0; i < iterMultiplier; i++)
(*iteratee) ();
difX = y[X_coord] - y[eqn0 + X_coord];
difY = y[Y_coord] - y[eqn0 + Y_coord];
dif = sqrt(difX * difX + difY * difY);
dif = sqrt(difX * difX + difY * difY);
ratio = dif / (numPixels * USTEP * diagonal);/* desired ratio = 1 */
if(printer >= 3)
PRINT "final ratio is now %lf; about 1 is target \n", ratio);
/*******
if(printer >= 3)
" J = %d:y1(x) = %11.11lf, y(x) = %11.11lf, dif=%11.11lf\n"
,J,y[eqn0+X_coord],y[X_coord],y[X_coord]-y[eqn0+X_coord]);
if(printer >= 2)
PRINT "diagonal = %lf final ratio = %lf \n",diagonal, ratio);
********/
/* now set ends of segment to be ya and yb */
for(i = zeroth; i < zeroth + vec_dim; i++)
ya[i] = yy[i];
for(i = zeroth; i < zeroth + vec_dim; i++)
yb[i] = y[i];
line_x1 = yy[X_coord];
line_x2 = y[X_coord];
line_y1 = yy[Y_coord];
line_y2 = y[Y_coord];
connect2 (y[eqn0 + X_coord], y[eqn0 + Y_coord], line_x1, line_y1);
/* connect2(line_x1,line_y1,line_x2,line_y2);*/
if(printer >= 2) {
"Now we compute images of ya-yb = (%lf,%lf)to(%lf,%lf) \n"
,line_x1, line_y1, line_x2, line_y2);
pause(3.);
}
}
LineForB(y1, y2)
double *y1,
*y2;
{
diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
+ (Y_upper - Y_lower) * (Y_upper - Y_lower));
line_x1 = y1[X_coord];
line_x2 = y2[X_coord];
line_y1 = y1[Y_coord];
line_y2 = y2[Y_coord];
stoplines(); /* for stopping drawing connected lines */
#ifdef X11
connectp(line_x1, line_y1);
#else
x_c = line_x1;
y_c = line_y1;
plot(x_c, y_c);
#endif
connectp(line_x2, line_y2);
return;
}
warning1 () { /* for lineForU */
{
"/n********************************** WARNING ******************************\n"
);
" THE FIRST STEP IS BIG. PERHAPS y1 IS NOT AT A PERIODIC ORBIT OF PERIOD NP\n"
);
" PROGRAM WILL ATTEMPT TO CARRY OUT COMMAND\n");
"********************************** WARNING ******************************\n\n"
);
}
}
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/