/* gcc -Ofast MandelJulia_X11.c -o mj_x11 -lX11 -lm */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>

#include <math.h>
#include <complex.h>

#include <X11/Xlib.h>
#include <X11/Xutil.h> // XPutPixel
#include <X11/cursorfont.h>


typedef uint32_t RGB32;


int MJ_ofsY = 100;
int M_ofsX = 100;
int J_ofsX = 860;


void image_Border(XImage *ximage) {
    RGB32 border = 0x555555;
    for (int x = 0; x < 640; x++) {
        XPutPixel(ximage, x,   0, border);
        XPutPixel(ximage, x, 479, border);
    }
    for (int y = 0; y < 480; y++) {
        XPutPixel(ximage,   0, y, border);
        XPutPixel(ximage, 639, y, border);
    }
}

void image_Mandelbrotset(XImage *ximage) {
    image_Border(ximage);

    for (int x = 1; x < 639; x++) {
        for (int y = 1; y < 479; y++) {
            int n = 0;
            complex double c = (x-420.0)/200.0 - I*(y-240.0)/200.0;
            complex double z = c;
            while (n < 1024 && cabs(z) <= 2.0) {
                z = z*z + c;
                n += 1;
            }
            RGB32 col = ((n * 4) & 0xFF) << 8;
            XPutPixel(ximage, x, y, col);
        }
    }
}

void image_Julia3(XImage *ximage, complex double c) {
    //image_Border(ximage);

    double bound = cabs(c);
    if (bound < 2.0) bound = 2.0;
    for (int x = 1; x < 639; x++) {
        for (int y = 1; y < 240; y++) { // 479
            int n = 0;
            complex double z = (x-320.0)/100.0 - I*(y-240.0)/100.0;
            while (n < 256 && cabs(z) <= bound) {
                z = z*z + c;
                n += 1;
            }
            RGB32 col = (n*8 & 0xFF) << 16;
            XPutPixel(ximage, x, y, col);
            XPutPixel(ximage, 639-x, 479-y, col);
        }
    }
}


//-------------------- external rays --------------------
// cf. https://abel.math.harvard.edu/~ctm/programs/
//     https://en.wikipedia.org/wiki/External_ray

complex double contsqrt(complex double z, complex double w) {
    complex double h = csqrt(z);
    if ( creal(h*conj(w)) < 0.0 ) {
        h = -h;
    }
    return h;
}

int len = 32;
int large = 200;
double res = 0.98;
double rayeps = 1e-4;
complex double riemann(double r, int p, int q, complex double z, complex double c) {
    complex double orb[len], norb[len];
    int n = 0;
    orb[0] = z;
    for (n = 0; n < len-1; n++) {
        if (cabs(orb[n]) > large) break;
        orb[n+1] = orb[n]*orb[n] + c;
        r *= 2.0;
        p = (2*p) % q;
    }
    norb[n] = pow(2.0, r) * cexp(I*2.0*M_PI*p/(double)q);
    while (n > 0) {
        norb[n-1] = contsqrt( norb[n] - c, orb[n-1] );
        n--;
    }
    return norb[0];
}

void draw_ray(XImage *ximage, int p, int q, complex double c) {
    double level = 2.0;
    double r = 8.0; // log_2(large)
    complex double z, z0;

    z = pow(2.0, r) * cexp(I*2.0*M_PI*p/(double)q);
    do {
        r *= res;
        z0 = z;
        z = riemann(r, p, q, z0, c);
        if (r < level*0.99) {
            int x = creal(z)*100;
            int y = cimag(z)*100;
            if (x > -320 && x < 319 && y > -239 && y < 240) {
                XPutPixel(ximage, 320+x, 240-y, 0xAAAA00);
            }
        }
    } while (r > rayeps);
}

void rational_rays(XImage *ximage, complex double c) {
    int q = 17;
    for (int p = 0; p < q; p++) {
        draw_ray(ximage, p, q, c);
    }
}

//
//-------------------- external rays --------------------

void image_Julia2(XImage *ximage, complex double c) {
    //image_Border(ximage);
    image_Julia3(ximage,  c);
    rational_rays(ximage, c);
}

void image_Julia1(XImage *ximage, complex double c) {
    //image_Border(ximage);
    for (int x = 1; x < 639; x++) {
        for (int y = 1; y < 479; y++) {
            XPutPixel(ximage, x, y, 0);
        }
    }

    int maxj = 1 << 14;
    int x = 0;
    int y = 0;

    complex double z = csqrt(0.25 - c) + 0.5;

    if (2.0*cabs(z) <= 1.0) {
        z = 1.0 - z;
    }

    for (int n = 0; n < maxj; n++) {

        z = csqrt( z - c );

        if ((rand() & 0xFF) < 128) {
            z = -z;
        }

        x = creal(z)*100;
        y = cimag(z)*100;

        if (x > -319 && x < 319 && y > -239 && y < 239)
        {
            RGB32 col = 0x0000AA;
            XPutPixel(ximage, 320+x, 240-y, col);
            XPutPixel(ximage, 320-x, 240+y, col);
        }
    }

    // critical point orbit
    int n = 0;
    complex double zc = 0;
    do {
        x = creal(zc) * 100;
        y = cimag(zc) * 100;
        if (x > -319 && x < 319 && y > -239 && y < 239) {
            XPutPixel(ximage, 320+x, 240-y, 0xAAAA00);
        }
        zc = zc*zc + c;
        n += 1;
    } while (n < 128 && cabs(zc) < 4);

}


int main (int argc, char *argv[]) {
    XEvent e;
    Atom wm_delete_window;

    Display *dsp = XOpenDisplay(NULL);

    Window w = XCreateSimpleWindow(dsp, DefaultRootWindow(dsp),
                                   100, 100, 1600, 800, 1,
                                   BlackPixel(dsp, DefaultScreen(dsp)),
                                   BlackPixel(dsp, DefaultScreen(dsp)) );

    XMapWindow(dsp, w);

    wm_delete_window = XInternAtom(dsp, "WM_DELETE_WINDOW", False);
    XSetWMProtocols(dsp, w, &wm_delete_window, 1);

    XSelectInput(dsp, w, StructureNotifyMask | ExposureMask | KeyPressMask | PointerMotionMask |
                         ButtonPressMask | Button1MotionMask);

    for ( e; e.type != MapNotify ; XWindowEvent(dsp, w, StructureNotifyMask, &e)) ;


    GC g = XCreateGC(dsp, w, 0, NULL);

    XSetForeground(dsp, g, 0x00FF00);

    XWindowAttributes wa = {0};
    XGetWindowAttributes(dsp, w, &wa);

    Cursor curs = XCreateFontCursor(dsp, XC_tcross);


    RGB32 mandel_px[640*480] = {0};
    XImage *img_M = XCreateImage(dsp, wa.visual, wa.depth, ZPixmap, 0, (char *)mandel_px, 640, 480, 32, 640*sizeof(RGB32));
    image_Mandelbrotset(img_M);
    XPutImage(dsp, w, g, img_M, 0, 0, M_ofsX, MJ_ofsY, 640, 480);

    RGB32 julia_px[640*480] = {0};
    XImage *img_J = XCreateImage(dsp, wa.visual, wa.depth, ZPixmap, 0, (char *)julia_px, 640, 480, 32, 640*sizeof(RGB32));
    image_Border(img_J);
    XPutImage(dsp, w, g, img_J, 0, 0, J_ofsX, MJ_ofsY, 640, 480);

    char emptyString[] = "                              ";
    char c_str[32] = "";
    RGB32 col = 0;

    int win_x = -1;
    int win_y = -1;
    uint32_t mask_return = 0;
    int jul = 0;

    int loop = 1;
    while (loop) {
        XNextEvent(dsp, &e);
        int _p = XPending(dsp);

        switch (e.type)
        {
            case Expose:
                    XPutImage(dsp, w, g, img_M, 0, 0, M_ofsX, MJ_ofsY, 640, 480);
                    XPutImage(dsp, w, g, img_J, 0, 0, J_ofsX, MJ_ofsY, 640, 480);
                    break;

            case ButtonPress:
            case MotionNotify:
                    win_x = e.xmotion.x;
                    win_y = e.xmotion.y;
                    mask_return = e.xmotion.state;
                    int mx = win_x - M_ofsX;
                    int my = win_y - MJ_ofsY;
                    if (mx > 0 && mx < 639 && my > 0 && my < 479) {
                        XDefineCursor(dsp, w, curs);

                        complex double c = (mx-420.0)/200.0 - I*(my-240.0)/200.0;

                        sprintf(c_str, "c = %+.3f%+.3fi", creal(c), cimag(c));
                        XSetForeground(dsp, g, 0x00AA00);
                        XDrawImageString(dsp, w, g, 200, 700, c_str, strlen(c_str));

                        if (e.xbutton.button == Button1 || (e.xmotion.state & Button1MotionMask)) jul = 1;
                        else
                        if (e.xbutton.button == Button3 || (e.xmotion.state & Button3MotionMask)) jul = 3;
                        else
                        if (e.xbutton.button == Button2 || (e.xmotion.state & Button2MotionMask)) jul = 2;
                        else jul = 0;

                        if ( jul && XPending(dsp) < 1 )
                        {
                            switch (jul) {
                                case 1:  image_Julia1(img_J, c); col = 0x0000CC; break;
                                case 2:  image_Julia2(img_J, c); col = 0xAAAA00; break;
                                case 3:  image_Julia3(img_J, c); col = 0xAA0000; break;
                            }
                            XPutImage(dsp, w, g, img_J, 0, 0, J_ofsX, MJ_ofsY, 640, 480);
                            XSetForeground(dsp, g, col);
                            XDrawImageString(dsp, w, g, 1000, 700, c_str, strlen(c_str));
                        }
                    }
                    else {
                        XDrawImageString(dsp, w, g, 200, 700, emptyString, 30);
                        XUndefineCursor(dsp, w);
                    }
                    break;

            case KeyPress:
                    loop = 0;
                    break;

            // close win: X connection to :0.0 broken (explicit kill or server shutdown).
            case ClientMessage:
                    if ((Atom)e.xclient.data.l[0] == wm_delete_window) {
                        loop = 0;
                    }
                    break;

            default:
                    XFlush(dsp);
                    break;
        }
    }

    XDestroyWindow(dsp, w);
    XCloseDisplay(dsp);

    return 0;
}