/* gcc -Wall hyperbole.c -lm -lX11 -lcairo */ #include #include #include #include #include #include #include #include #include #include #include int fps = 20; #define WIDTH (512) #define HEIGHT (512) cairo_surface_t *cs; cairo_surface_t *csout; void paint(int curtime) { cairo_t *c; int i, j; c=cairo_create(cs); /* clear screen */ cairo_rectangle(c, 0, 0, WIDTH, HEIGHT); cairo_set_source_rgb(c, 0.8, 0.8, 0.7); cairo_fill(c); cairo_set_source_rgb(c, 0.2, 0.5, 0.7); for (i = 0; i < 100; i++) for (j = 0; j < 100; j++) { double m[2][2]; double x = (i-50) / 10.; double y = (j-50) / 10.; double t = curtime; #if 0 /* rotation */ double theta = 2 * M_PI * t/100; m[0][0] = cos(theta); m[0][1] = -sin(theta); m[1][0] = sin(theta); m[1][1] = cos(theta); #endif #if 0 /* exp */ t = (int)t % 100; m[0][0] = exp(t/100); m[0][1] = 0; m[1][0] = 0; m[1][1] = exp(-t/100); double theta = M_PI/4 + M_PI/4/3; m[0][0] = exp(t/100) * cos(theta); m[0][1] = exp(t/100) * -sin(theta); m[1][0] = exp(-t/100) * sin(theta); m[1][1] = exp(-t/100) * cos(theta); #endif #if 1 /* exp */ t = (int)t % 100; m[0][0] = exp(t/100); m[0][1] = 0; m[1][0] = 0; m[1][1] = exp(-t/100); double theta = M_PI/4 + M_PI/4/3; double costheta = (sqrt(5)-1) / 2; double sintheta = sqrt(1 - costheta * costheta); m[0][0] = exp(t/100) * costheta; m[0][1] = exp(t/100) * -sintheta; m[1][0] = exp(-t/100) * sintheta; m[1][1] = exp(-t/100) * costheta; #endif #if 0 /* t 1/t */ t = (int)t % 100; m[0][0] = (t+1)/100; m[0][1] = 0; m[1][0] = 0; m[1][1] = 1/m[0][0]; double theta = M_PI/4 + M_PI/4/3; m[0][0] = exp(t/100) * cos(theta); m[0][1] = exp(t/100) * -sin(theta); m[1][0] = exp(-t/100) * sin(theta); m[1][1] = exp(-t/100) * cos(theta); #endif #if 0 /* m = (a b * c d) * * x' = ax + by * y' = cx + dy */ /* linear interpolation - (1,0) goes to (1,1) and (1,1) to (2,3) */ t = (int)t % 100; m[0][0] = 1; m[0][1] = 1 * t/100; m[1][0] = 1 * t/100; m[1][1] = 1 + t/100; #endif #if 0 /* * 1 = a * 0 = a + b * b = -1 * 1 = c * 1 = c + d * d = 0 * so (1 0 goes to (1 -1 * 0 1) 1 0) */ /* linear interpolation - (1,0) goes to (1,1) and (1,1) to (0,1) */ t = (int)t % 100; m[0][0] = 1; m[0][1] = -t/100; m[1][0] = t/100; m[1][1] = 1 - t/100; #endif #if 0 /* sinh(x) = (e^x - e^-x) / 2 * cosh(x) = (e^x + e^-x) / 2 * * * * * * * * * * * * D C * * * A B * * * * * * * * * * C * * * D * * * * B * * * A * * * * * * * (1,0) goes to (1,1) and (1,1) to (2,3) with unit surface (1,1)->(2,3) linear interpolation */ t = (int)t % 200; double ax = 0; double ay = 0; double bx = 1; double by = (exp(t/100)-exp(-t/100))/2; double cx = 1 + t/100; double cy = 1 + 2*t/100; bx = (exp(t/100) + exp(-t/100)) / 2; by = (exp(t/100) - exp(-t/100)) / 2; cx = (exp((t+100)/100) + exp(-(t+100)/100)) / 2; cy = (exp((t+100)/100) - exp(-(t+100)/100)) / 2; double px = (ax*(by-ay)*(by-ay) + cx*(bx-ax)*(bx-ax) -(bx-ax)*(by-ay)*(ay-cy)) / ((by-ay)*(by-ay) + (bx-ax)*(bx-ax)); double py = ay - (ay-by)*(px-ax) / (bx-ax); double k = ((px-cx)*(px-cx)+(py-cy)*(py-cy)) * ((bx-ax)*(bx-ax)+(by-ay)*(by-ay));; k = 1/k; double bbx = ax + k * (bx-ax); double bby = ay + k * (by-ay); bx = bbx; by = bby; /* so (1,0) to (bx,by) and (1,1) to (cx,cy) * x' = ax + by * y' = cx + dy * * bx = a * cx = a + b * by = c * cy = c + d */ { double a = bx; double b = cx - a; double c = by; double d = cy - c; m[0][0] = a; m[0][1] = b; m[1][0] = c; m[1][1] = d; } #endif #if 1 /* https://observablehq.com/@chnn/the-modular-flow-in-tex-mathcal-l-and-tex-mathbb-h#latticeInOrbit */ /* http://www.josleys.com/articles/ams_article/Lorenz3.htm#orbits */ /* A : 1 1 1 2 find the eigenvalues l (det(A-lI) = 0): A - lI : 1-l 1 1 2-l det(A-lI) = (1-l)(2-l) - 1 = 0 l = (3 +/- sqrt(5)) / 2 2 eigenvalues: l1 = (3-sqrt(5)) / 2 l2 = (3+sqrt(5)) / 2 find eigenvectors: (A - lI)x = 0 x = (x1 x2) (1-l)x1 + x2 = 0 x1 + (2-l)x2 = 0 gaussian elimiation(?) (this I don't remember very well, (1-l)(2-l) = 1 because det(A-lI) = 0) (1-l)x1 + x2 = 0 (1-l)x1 + (1-l)(2-l) x2 = 0 so we get: x2 = (l-1)x1 we set x1 = 1 and get x2 for l1 and l2 then "glue" eigenvectors to get the matrix U: U: 1 1 l1-1 l2-1 find U-1 UU-1 = I a b c d 1 1 a+c b+d l1-1 l2-1 (l1-1)a + (l2-1)c (l1-1)b + (l2-1)d a+c = 1 b+d = 0 (l1-1)a + (l2-1)c = 0 (l1-1)b + (l2-1)d = 1 a=1-c (l1-1)(1-c) + (l2-1)c = 0 (l1-1) + c (l2-l1) = 0 c = (1-l1) / (l2-l1) b = -d -(l1-1)d + (l2-1)d = 1 d = 1 / (l2-l1) U-1 : 1 - (1-l1) / (l2-l1) -1/(l2-l1) (1-l1) / (l2-l1) 1 / (l2-l1) that is: (l2-1) / (l2-l1) -1/(l2-l1) (1-l1) / (l2-l1) 1 / (l2-l1) we get rid of (l2-l1) in code below (1 0) goes to (l2-1 -(l1-1)) (1.618034 0.618034) (0 1) goes to (-1 1) we need to rotate for dots to be vertically/horizontally aligned. For that, we visually find the dots that map to A and B as seen in the following graph (minus rotation, which we want to get rid of) for the right t (to be found too): * * * * * * * A * * * * O B * * * * * * * * * * * We find that it's for (1, 1) and (0, -1) Using the transform, we have: 1 1 => (l2-2) e(t) (2-l1) e(-t) 0 -1 => 1 e(t) -1 e(-t) We want to find t for which it's orthognal. We want orthogonality, so scalar product is = 0: (l2-2)e(2t) + (l1-2)e(-2t) = 0 (l2-2)e(4t) + (l1-2) = 0 e(4t) = (2-l1)/(l2-2) We find t: t = ln((2-l1)/(l2-2)) / 4 Using (0, -1) we can find the cos of the angle (_cos in the code below): (it's x / sqrt(x*x+y*y) where (x, y) is the destination point of (0,-1)) (0 -1) at t = ln(...): e(ln((2-l1)/(l2-2)) / 4) e(-ln((2-l1)/(l2-2)) / 4) Then we get m2, inverse rotation matrix: _cos -_sin _sin _cos */ #endif double l1 = (3-sqrt(5)) / 2; double l2 = (3+sqrt(5)) / 2; if (x == 0 && y == -.1) { printf("%g %g\n", exp(1./4) * (2-l1)/(l2-2), -exp(-1./4) * (2-l1)/(l2-2)); printf("%g %g\n", exp(log((2-l1)/(l2-2)) / 4), -exp(-log((2-l1)/(l2-2)) / 4)); printf("cos(x)=%g\n", exp(log((2-l1)/(l2-2)) / 4) / sqrt(exp(log((2-l1)/(l2-2)) / 4)*exp(log((2-l1)/(l2-2)) / 4) + exp(-log((2-l1)/(l2-2)) / 4)*exp(-log((2-l1)/(l2-2)) / 4))); } double m2[2][2]; double _cos = exp(log((2-l1)/(l2-2)) / 4) / sqrt(exp(log((2-l1)/(l2-2)) / 4)*exp(log((2-l1)/(l2-2)) / 4) + exp(-log((2-l1)/(l2-2)) / 4)*exp(-log((2-l1)/(l2-2)) / 4)); double _sin = sqrt(1 -_cos*_cos); m2[0][0] = _cos; m2[0][1] = -_sin; m2[1][0] = _sin; m2[1][1] = _cos; /* loop period is log(l1) or log(l2), so scale t from [0..1] to [0..log(l2)] */ t = (int)t % 100; double et1 = exp(t/100 * log(l2)); double et2 = exp(-t/100 * log(l2)); t = (int)t % 30; /* start at right t [log((2-l1)/(l2-2)) / 4] to have horizontal/vertical at first image */ t = t/30*log(l1) + log((2-l1)/(l2-2)) / 4; //t = log((2-l1)/(l2-2)) / 4; et1 = exp(t); et2 = exp(-t); /* the magic is there - multiply U-1 by (et1, et2), which may also scale things * (knowing that we got rid of "/(l2-l1)") * anyway I don't get it fully clear in my head * it works, that's all I can say */ m[0][0] = et1 * (l2-1); m[0][1] = et1 * -1; m[1][0] = et2 * (1-l1); m[1][1] = et2 * 1; double x3 = m[0][0] * x + m[0][1] * y; double y3 = m[1][0] * x + m[1][1] * y; /* inverse rotate */ double x2 = m2[0][0] * x3 + m2[0][1] * y3; double y2 = m2[1][0] * x3 + m2[1][1] * y3; //x2 /= (l2-l1); //y2 /= (l2-l1); if (x == 0 && y == 0) printf("[%g %g] [%g %g] %g\n", m[0][0], m[0][1], m[1][0], m[1][1], m[0][0] * m[0][1] + m[1][0] * m[1][1]); if (x == 0 && y == 0 || x == .1 && y == 0 || x == .0 && y == .1 || x == .0 && y == -.1) printf("%g %g -> %g %g\n", x, y, x2, y2); if (x == 0 && y == 0 || x == .1 && y == 0 || x == .0 && y == .1) cairo_set_source_rgb(c, 0.8, 0.4, 0.5); else cairo_set_source_rgb(c, 0.2, 0.5, 0.7); if(x == .0 && y == .1) cairo_set_source_rgb(c, 0.3, 0.9, 0.5); if(x == .1 && y == .1) cairo_set_source_rgb(c, 0.3, 0.2, 0.5); if(x == 0 && y == -.1) cairo_set_source_rgb(c, 0., 0., 0.9); static double px1=0, py1=0, px2=0, py2=0; if(x == .1 && y == .1) { px1 = x2; py1 = y2; } if(x == 0 && y == -.1) { px2 = x2; py2 = y2; } if(x==1&&y==1)printf("%g\n", px1*px2+py1*py2); x = x2*256+256; y = 511-(y2*256+256); // cairo_rectangle(c, x, y, 10, 10); cairo_arc(c, x, y, 16, 0, 2*M_PI); cairo_fill(c); } cairo_set_source_rgba(c, 0.1, 0.1, 0.1, .4); if (0) for (i = 0; i < 100; i++) for (j = 0; j < 100; j++) { double m[2][2]; double x = (i-50) / 10.; double y = (j-50) / 10.; x = x*256+256; y = 511-(y*256+256); cairo_rectangle(c, x, y, 10, 10); cairo_fill(c); } cairo_destroy(c); /* and put it all to the xpixmap */ c=cairo_create(csout); cairo_set_source_surface(c, cs, 0, 0); cairo_rectangle(c, 0, 0, WIDTH, HEIGHT); cairo_paint(c); cairo_destroy(c); //getchar(); } void expose(Display *d, Window w, Pixmap p) { XCopyArea(d, p, w, DefaultGC(d, DefaultScreen(d)), 0,0, WIDTH, HEIGHT, 0,0); } int main(int n, char **v) { Display *d; Window w; Pixmap p; fd_set r; int repaint; int redraw; int i; int ret; struct timeval time_out; int curtime = 0; for (i = 1; i < n; i++) { } cs = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, WIDTH, HEIGHT); /* x window mode */ d = XOpenDisplay(0); if (!d) abort(); w = XCreateSimpleWindow(d, DefaultRootWindow(d), 0, 0, WIDTH, HEIGHT, 0, 0, 0); XSelectInput(d, w, KeyPressMask | ExposureMask); XMapWindow(d, w); p = XCreatePixmap(d, w, WIDTH,HEIGHT, DefaultDepth(d, DefaultScreen(d))); csout = cairo_xlib_surface_create(d, p, DefaultVisual(d, DefaultScreen(d)), WIDTH, HEIGHT); repaint = 1; redraw = 0; while (1) { XEvent ev; while (XPending(d)) { XNextEvent(d, &ev); switch (ev.type) { case Expose: redraw = 1; break; case KeyPress: { char buf[65]; KeySym ks; XLookupString(&ev.xkey, buf, 64, &ks, NULL); switch (ks) { case XK_Right: case 'n': printf("n"); fflush(stdout); break; case XK_Left: case 'p': printf("p"); fflush(stdout); break; case XK_Escape: case 'q': exit(1); } break; } default: break; } } if (repaint) { paint(curtime); redraw = 1; repaint = 0; } if (redraw) { expose(d, w, p); redraw = 0; } FD_ZERO(&r); FD_SET(ConnectionNumber(d), &r); XFlush(d); /* this time keeping is wrong, we sleep too much because we don't * take into account computation time of previous frame. */ time_out.tv_sec = 0; time_out.tv_usec = 1000000/fps; ret = select(ConnectionNumber(d)+1, &r, 0, 0, &time_out); if (ret == -1) abort(); if (ret == 0) { repaint = 1; curtime++; } } return 0; }