Physics 版 (精华区)
发信人: zjliu (秋天的萝卜), 信区: Physics
标 题: (转载) 混沌吸引子的作图程序源码
发信站: 哈工大紫丁香 (2003年05月21日22:17:35 星期三), 站内信件
【 以下文字转载自 Lixueyuan 讨论区 】
【 原文由 zjliu 所发表 】
发信站: BBS 水木清华站
下面这个程序是从一个Hacker站点得来的.适用于X-window系统.我已经运行通过.
画出logistic映射的混沌情况.
/*#if 0
cc -o chaos chaos.c -lm -lX11
exit
#endif
*/
/*
* Chaos program. Copyright 1991 Ken Shirriff
* Send comments, etc. to shirriff@sprite.Berkeley.EDU
*
* This program draws the "chaotic attractor".
* Usage: chaos [-r min max] [-i ilow ihigh] [-s step]\n");
* -r gives the minimum and maximum values of the logistic equation parameter.
* -i gives the iteration range to be plotted and used for the Lyapunov
* exponent computation. (The iterations below ilow are not used, so they
* exponent computation. (The iterations below ilow are not used, so they
* can be used to allow the function to settle.
* -s gives the step size (i.e. how big the steps are in logistic equation
* parameter space. 1=1 pixel.
*
* After drawing, left click prints the X position of the mouse in terms
* of the logistic equation parameter.
* Right click ends the program.
*
* Compile program with "sh chaos.c"
*
*/
#include <X11/Xlib.h>
#include <stdio.h>
#include <math.h>
struct Point {
double x,y;
};
typedef struct Point pos_t;
/*
* SIZE is the size of the X window to use.
*/
#define SIZE 600
/*
* The logistic equation.
*/
double logis(x,s)
double x,s;
{
double xp;
xp = s*x*(1-x);
return xp;
}
main(argc,argv)
int argc;
char **argv;
{
double s; /* Our logistics equation driving parameter. */
int i;
int i;
double x; /* Our iterated x value. */
pos_t p;
float lya; /* Lyapunov exponent. */
int N0=10, NT=100; /* Low and high iterations to plot. */
float S0=0, S1=4.0; /* Low and high parameter values. */
float div=1; /* Plotting resolution in pixels. */
#define abs(x) ((x)>0?(x):-(x))
if (argc < 1) {
usage:
fprintf(stderr,"Usage: [-r min max] [-i ilow ihigh] [-s step]\n");
exit(1);
}
while (1) {
if (argc==1) {
break;
} else if (strncmp("-h",argv[1],2)==0) {
goto usage;
} else if (strcmp("-r",argv[1])==0 && argc>3) {
sscanf(argv[2],"%f", &S0);
sscanf(argv[3],"%f", &S1);
if (S0>=S1) {
if (S0>=S1) {
fprintf(stderr,"Invalid range\n");
exit(-1);
}
argc -= 3;
argv += 3;
} else if (strcmp("-i",argv[1])==0 && argc>3) {
sscanf(argv[2],"%d", &N0);
sscanf(argv[3],"%d", &NT);
argc -= 3;
argv += 3;
} else if (strcmp("-s",argv[1])==0 && argc>2) {
sscanf(argv[2],"%f", &div);
argc -= 2;
argv += 2;
} else if (strcmp("-i",argv[1])==0 && argc>3) {
sscanf(argv[2],"%d", &N0);
sscanf(argv[3],"%d", &NT);
argc -= 3;
argv += 3;
} else {
goto usage;
}
}
}
xinit();
range(S0,-2.,S1,1.5);
/*
* Loop across the different 'C' values (in the variable s)
*/
for (s=S0;s<S1;s+=(S1-S0)/SIZE*div) {
x = .5;
lya = 0;
/*
* Do NT iterations of the logistics equation.
*/
for (i=0;i<NT;i++) {
x = logis(x,s);
if (i>N0) {
/*
* Plot the point.
*/
p.x = s;
p.y = x;
p.y = x;
point(&p);
/*
* Sum up the Lyapunov exponent.
*/
lya += log(abs(s-2*x*s))/(NT-N0);
}
}
/*
* Plot the Lyapunov exponent and 0 axis.
*/
p.x = s;
p.y = lya/2-1; /* Scale lya to lower part of window. */
point(&p);
p.x = s;
p.y = 0./2-1;
point(&p);
}
/*
* Flush the drawing buffer.
*/
drawflush();
drawflush();
printf("Done\n");
click();
endit();
}
/*
* -------------------------------------------------------------
* X drivers
*/
Display *dp;
Window w;
GC gc;
float Xxmin, Xxmax, Xymin, Xymax;
#define SCALEX(v) ((int)(((v)-(Xxmin))*SIZE/(Xxmax-(Xxmin))+.5))
#define ISCALEX(n) (Xxmin+(n)*(Xxmax-Xxmin)/SIZE)
#define SCALEY(v) (SIZE-((int)(((v)-(Xymin))*SIZE/(Xymax-(Xymin))+.5)))
#define ISCALEY(n) (Xymin+(SIZE-n)*(Xymax-Xymin)/SIZE)
/*
* Initialize the X window.
*/
xinit()
{
Window root;
Screen *sc;
int dscreen;
XGCValues gcvals;
XSetWindowAttributes watt;
int depth;
XEvent xevent;
int done=0;
if ((dp=XOpenDisplay(NULL))==NULL) {
fprintf(stderr,"Could not open Display!0\n");
exit(1);
}
XSynchronize(dp, True);
dscreen = XDefaultScreen(dp);
sc = ScreenOfDisplay(dp,dscreen);
sc = ScreenOfDisplay(dp,dscreen);
depth = DefaultDepth(dp,dscreen);
root = DefaultRootWindow(dp);
watt.background_pixel = WhitePixelOfScreen(sc);
watt.bit_gravity = StaticGravity;
w = XCreateWindow(dp,root,0,0,SIZE,SIZE,0,
depth,InputOutput,CopyFromParent,CWBackPixel|CWBitGravity, &watt);
XStoreName(dp,w,"Test");
gcvals.background = BlackPixelOfScreen(sc);
gcvals.foreground = BlackPixelOfScreen(sc);
gc = XCreateGC(dp,w, GCForeground|GCBackground, &gcvals);
XSetFunction(dp,gc,GXcopy);
XSelectInput(dp,w,ExposureMask|EnterWindowMask|PointerMotionMask|
ButtonPressMask|StructureNotifyMask|SubstructureNotifyMask);
XMapWindow(dp,w);
XSync(dp,False);
XClearArea(dp,w,0,0,SIZE, SIZE, False);
while (!done) {
XNextEvent(dp, &xevent);
switch (((XAnyEvent *)&xevent)->type) {
case ButtonPress:
if ( ((XButtonEvent*)&xevent)->button == Button3) {
done = 1;
done = 1;
}
break;
case Expose:
done = 1;
break;
default:
break;
}
}
}
/*
* Set the plotting range.
*/
range(pxmin, pymin, pxmax, pymax)
float pxmin, pxmax, pymin, pymax;
{
Xxmin = pxmin;
Xxmax = pxmax;
Xymin = pymin;
Xymax = pymax;
}
下面还有喔 (76%) │ 结束 ← <q> │ ↑/↓/PgUp/PgDn 移动 │ ? 辅助说明 │
}
/*
* Plot the point. Points are buffered and plotted in large bunches
* to make it run faster.
*/
#define BUF 1000
XPoint buf[BUF];
int bufp=0;
point(p)
pos_t *p;
{
if (bufp>=BUF) {
XDrawPoints(dp,w,gc,buf,bufp,CoordModeOrigin);
bufp=0;
}
buf[bufp].x = SCALEX(p->x);
buf[bufp].y = SCALEY(p->y);
bufp++;
}
/*
* Flush the point drawing buffer.
* Flush the point drawing buffer.
*/
drawflush()
{
if (bufp>0) {
XDrawPoints(dp,w,gc,buf,bufp,CoordModeOrigin);
bufp=0;
}
}
/*
* Wait for mouse click.
*/
click()
{
XEvent xevent;
int done=0;
while (!done) {
XNextEvent(dp, &xevent);
switch (((XAnyEvent *)&xevent)->type) {
case ButtonPress:
if ( ((XButtonEvent*)&xevent)->button == Button3) {
done = 1;
done = 1;
} else if ( ((XButtonEvent*)&xevent)->button == Button1) {
/*
* If we get a click, print the mouse position.
*/
int x,y;
float fy;
x = ((XButtonEvent *)&xevent)->x;
y = ((XButtonEvent *)&xevent)->y;
fy = ISCALEY(y);
/*
* Rescale the Lyapunov window coordinates.
*/
if (fy<-.2) {
printf("fy=%f\n", fy);
fy = (fy+1)*2;
}
printf("(%f,%f)\n", ISCALEX(x), fy);
}
break;
default:
break;
}
printf("fy=%f\n", fy);
fy = (fy+1)*2;
}
printf("(%f,%f)\n", ISCALEX(x), fy);
}
break;
default:
break;
} }
}
/*
* Close the window.
*/
endit()
{
XDestroyWindow(dp, w);
XCloseDisplay(dp);
}
--
※ 来源:．哈工大紫丁香 http://bbs.hit.edu.cn [FROM: 202.118.229.86]
--
※ 转载:．哈工大紫丁香 bbs.hit.edu.cn．[FROM: 202.118.229.86]
Powered by KBS BBS 2.0 (http://dev.kcn.cn)
页面执行时间：202.436毫秒