//===### The file — 'ContoArrayJump15_test.cc'===
//<pre>
// ContoArrayJump15_test.cc
// Self-contained test program with ContoArrayJump15 contour routine.
// Produces identity_map.eps, sqrt_map.eps, log_map.eps
// Compile: g++ -O2 -std=gnu++98 ContoArrayJump15_test.cc -o ContoArrayJump15_test

#include <cstdio>
#include <cmath>
#include <complex>
#include <cstdlib>

typedef std::complex<double> cd;

#ifndef EPS_ZERO
#define EPS_ZERO 1e-12
#endif

// --------------------- ado (EPS prolog helper) ---------------------
void ado(FILE* out, double W, double H) {
    if (!out) return;
    std::fprintf(out, "%%!PS-Adobe-3.0 EPSF-3.0\n");
    std::fprintf(out, "%%%%BoundingBox: 0 0 %.0f %.0f\n", W, H);
    std::fprintf(out, "/M { moveto } def\n");
    std::fprintf(out, "/L { lineto } def\n");
    std::fprintf(out, "/W { setlinewidth } def\n");
    std::fprintf(out, "/RGB { setrgbcolor } def\n");
    std::fprintf(out, "/S { stroke } def\n");
    std::fprintf(out, "1 setlinejoin 1 setlinecap\n");
}

// ------------------ ContoArrayJump15 (flattened array) ------------------
// field: flattened array of size (M+1)*(N+1), index = i + j*(M+1)
// X[0..M], Y[0..N]
// max_jump: skip edge if fabs(fa-fb) > max_jump; set <0 to disable
void ContoArrayJump15(FILE* out,
                      double *field, double *X, double *Y,
                      int M, int N,
                      double level,
                      double max_jump)
{
    // M = number of cells horizontally; NX = nodes = M+1
    int NX = M + 1;
    int NY = N + 1;

    for (int j = 0; j < N; ++j) {
        for (int i = 0; i < M; ++i) {

            // corners: bottom-left (f00), bottom-right (f10),
            //          top-right (f11), top-left (f01)
            double f00 = field[i + j*NX];
            double f10 = field[(i+1) + j*NX];
            double f11 = field[(i+1) + (j+1)*NX];
            double f01 = field[i + (j+1)*NX];

            // skip cells that contain NaN or huge values
            if (!std::isfinite(f00) || !std::isfinite(f10) || !std::isfinite(f11) || !std::isfinite(f01))
                continue;

            double vmin = f00, vmax = f00;
            if (f10 < vmin) vmin = f10;
            if (f11 < vmin) vmin = f11;
            if (f01 < vmin) vmin = f01;
            if (f10 > vmax) vmax = f10;
            if (f11 > vmax) vmax = f11;
            if (f01 > vmax) vmax = f01;

            // quick reject
            if (level < vmin - EPS_ZERO || level > vmax + EPS_ZERO) continue;

            // Per-edge tests and intersections
            // edge order we'll record: 0=bottom,1=right,2=top,3=left
            double xp[4], yp[4];
            bool has[4] = {false,false,false,false};
            int cnt = 0;

            // bottom edge: (x0,y0) -> (x1,y0) with values f00,f10
            {
                double fa = f00, fb = f10;
                bool signchange = (fa - level)*(fb - level) <= 0.0;
                bool jumpok = (max_jump < 0.0) || (fabs(fa - fb) <= max_jump);
                if (signchange && jumpok) {
                    double x0 = X[i], x1 = X[i+1], y0 = Y[j];
                    double denom = (fb - fa);
                    double t = (fabs(denom) < EPS_ZERO) ? 0.5 : ((level - fa) / denom);
                    if (t < -1e6) t = 0.5; // safety
                    xp[cnt] = x0 + t*(x1 - x0);
                    yp[cnt] = y0;
                    has[cnt] = true;
                    ++cnt;
                }
            }

            // right edge: (x1,y0) -> (x1,y1) with values f10,f11
            {
                double fa = f10, fb = f11;
                bool signchange = (fa - level)*(fb - level) <= 0.0;
                bool jumpok = (max_jump < 0.0) || (fabs(fa - fb) <= max_jump);
                if (signchange && jumpok) {
                    double x0 = X[i+1], y0 = Y[j], y1 = Y[j+1];
                    double denom = (fb - fa);
                    double t = (fabs(denom) < EPS_ZERO) ? 0.5 : ((level - fa) / denom);
                    xp[cnt] = x0;
                    yp[cnt] = y0 + t*(y1 - y0);
                    has[cnt] = true;
                    ++cnt;
                }
            }

            // top edge: (x0,y1) -> (x1,y1) with values f01,f11
            {
                double fa = f01, fb = f11;
                bool signchange = (fa - level)*(fb - level) <= 0.0;
                bool jumpok = (max_jump < 0.0) || (fabs(fa - fb) <= max_jump);
                if (signchange && jumpok) {
                    double x0 = X[i], x1 = X[i+1], y0 = Y[j+1];
                    double denom = (fb - fa);
                    double t = (fabs(denom) < EPS_ZERO) ? 0.5 : ((level - fa) / denom);
                    xp[cnt] = x0 + t*(x1 - x0);
                    yp[cnt] = y0;
                    has[cnt] = true;
                    ++cnt;
                }
            }

            // left edge: (x0,y0) -> (x0,y1) with values f00,f01
            {
                double fa = f00, fb = f01;
                bool signchange = (fa - level)*(fb - level) <= 0.0;
                bool jumpok = (max_jump < 0.0) || (fabs(fa - fb) <= max_jump);
                if (signchange && jumpok) {
                    double x0 = X[i], y0 = Y[j], y1 = Y[j+1];
                    double denom = (fb - fa);
                    double t = (fabs(denom) < EPS_ZERO) ? 0.5 : ((level - fa) / denom);
                    xp[cnt] = x0;
                    yp[cnt] = y0 + t*(y1 - y0);
                    has[cnt] = true;
                    ++cnt;
                }
            }

            // Now analyze intersections:
            if (cnt == 2) {
                // Simple case: draw between the two points
                std::fprintf(out,"newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n",
                             xp[0], yp[0], xp[1], yp[1]);
                continue;
            }
            else if (cnt == 0) {
                // no edge crossings due to jumps or exact corner cases: handle corner equalities
                // Handle special corner-equality cases using EPS_ZERO
                bool c00 = fabs(f00 - level) < EPS_ZERO;
                bool c10 = fabs(f10 - level) < EPS_ZERO;
                bool c11 = fabs(f11 - level) < EPS_ZERO;
                bool c01 = fabs(f01 - level) < EPS_ZERO;
                int eqcount = (c00?1:0)+(c10?1:0)+(c11?1:0)+(c01?1:0);
                if (eqcount >= 2) {
                    // draw a small diagonal to indicate the level passes near corners
                    double xa = X[i], xb = X[i+1], ya = Y[j], yb = Y[j+1];
                    // choose a reasonable diagonal if opposite corners equal
                    if ((c00 && c11) || (c10 && c01)) {
                        std::fprintf(out,"newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n", xa, ya, xb, yb);
                    } else if (c00 || c11) {
                        std::fprintf(out,"newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n", xa, ya, xb, yb);
                    } else {
                        std::fprintf(out,"newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n", xa, yb, xb, ya);
                    }
                }
                continue;
            }
            else if (cnt == 4) {
                // Ambiguous saddle case: use bilinear-center decider
                // compute center value (cheap but effective)
                double fC = 0.25*(f00 + f10 + f11 + f01);
                bool centerHigh = (fC >= level);

                // decide diagonals:
                // diagonal1: bottom <-> top (xp[0],yp[0]) <-> (xp[2],yp[2]) [bottom & top]
                // diagonal2: left <-> right (xp[3],yp[3]) <-> (xp[1],yp[1]) [left & right]
                // But we need to be careful mapping xp positions to edges:
                // Under our filling order, xp[0]=bottom, xp[1]=right, xp[2]=top, xp[3]=left if all present.
                // Determine which diagonal connects the two 'high' corners.
                // Build corner boolean for >=level (with EPS):
                bool b00 = (f00 >= level - EPS_ZERO);
                bool b10 = (f10 >= level - EPS_ZERO);
                bool b11 = (f11 >= level - EPS_ZERO);
                bool b01 = (f01 >= level - EPS_ZERO);
                // If high corners are (00 and 11) -> connect bottom-top diagonal (xp[0] <-> xp[2]).
                // If high corners are (10 and 01) -> connect left-right diagonal (xp[3] <-> xp[1]).
                bool high00_11 = (b00 && b11 && !b10 && !b01);
                bool high10_01 = (b10 && b01 && !b00 && !b11);

                if (high00_11 || (!high10_01 && centerHigh)) {
                    // connect bottom <-> top
                    std::fprintf(out,"newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n", xp[0], yp[0], xp[2], yp[2]);
                } else {
                    // connect left <-> right
                    std::fprintf(out,"newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n", xp[3], yp[3], xp[1], yp[1]);
                }
                continue;
            }
            else {
                // cnt == 3 or other degenerate combinations: connect pairs closest together
                // Find two points with maximal separation (roughly drawing the main chord)
                double bestd = -1.0;
                int a=0,b=1;
                for (int p=0;p<cnt;p++) for(int q=p+1;q<cnt;q++){
                    double dx = xp[p]-xp[q], dy = yp[p]-yp[q];
                    double d2 = dx*dx+dy*dy;
                    if (d2 > bestd) { bestd = d2; a=p; b=q; }
                }
                std::fprintf(out,"newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n", xp[a], yp[a], xp[b], yp[b]);
            }
        }
    }
}

// -------------------------- helpers & main tests --------------------------
void draw_grid_eps(FILE* out, double *X, double *Y, int M, int N) {
    // draw light grid lines between integer ticks
    std::fprintf(out, "0 0 0 RGB .006 W\n");
    // vertical lines at integer x (if inside)
    double xmin = X[0], xmax = X[M];
    double ymin = Y[0], ymax = Y[N];
    for (int xi = (int)ceil(xmin); xi <= (int)floor(xmax); ++xi) {
        std::fprintf(out,"%.12g %.12g M %.12g %.12g L S\n", (double)xi, ymin, (double)xi, ymax);
    }
    for (int yj = (int)ceil(ymin); yj <= (int)floor(ymax); ++yj) {
        std::fprintf(out,"%.12g %.12g M %.12g %.12g L S\n", xmin, (double)yj, xmax, (double)yj);
    }
}

int main(int argc, char** argv) {
    // grid parameters (cells)
    const int M = 200;   // number of cells horizontally
    const int N = 200;   // number of cells vertically
    const double xmin = -3.0, xmax = 3.0;
    const double ymin = -3.0, ymax = 3.0;
    const double dx = (xmax - xmin) / double(M);
    const double dy = (ymax - ymin) / double(N);
    const int NX = M+1, NY = N+1;

    // allocate grid coords
    double *X = (double*)malloc(sizeof(double)*NX);
    double *Y = (double*)malloc(sizeof(double)*NY);

    for (int i=0;i<NX;i++) X[i] = xmin + i*dx;
    for (int j=0;j<NY;j++) Y[j] = ymin + j*dy;

    // allocate fields (flattened)
    double *field = (double*)malloc(sizeof(double)*NX*NY);

    // ---------- 1) identity map (Re and Im)
    // Re map
    for (int j=0;j<NY;j++) for (int i=0;i<NX;i++) field[i + j*NX] = X[i];
    FILE *o = fopen("identity_map.eps","w");
    ado(o, 422.0, 422.0);
    std::fprintf(o,"211 211 translate\n 100 100 scale\n");
    draw_grid_eps(o, X, Y, M, N);
    // draw Im=integer levels as horizontal colored lines: negative red, zero magenta, positive blue
    for (int k=-3;k<=3;k++) {
        if (k < 0) std::fprintf(o,"1 0 0 RGB .03 W\n"); else if (k>0) std::fprintf(o,"0 0 1 RGB .03 W\n"); else std::fprintf(o,".6 0 .6 RGB .04 W\n");
        ContoArrayJump15(o, field, X, Y, M, N, (double)k, -1.0); // max_jump<0 => disabled
    }
    // Re levels vertical thin
    for (int k=-3;k<=3;k++) {
        if (k < 0) std::fprintf(o,"1 0 0 RGB .01 W\n"); else if (k>0) std::fprintf(o,"0 0 1 RGB .01 W\n"); else std::fprintf(o,".5 0 .5 RGB .02 W\n");
        ContoArrayJump15(o, field, X, Y, M, N, (double)k, -1.0);
    }
    std::fprintf(o,"showpage\n%%%%Trailer\n");
    fclose(o);

    // ---------- 2) sqrt map (principal branch) -> plot Im and Re contours
    for (int j=0;j<NY;j++) for (int i=0;i<NX;i++) {
        cd z(X[i], Y[j]);
        cd w = std::sqrt(z); // principal
        field[i + j*NX] = std::real(w); // we'll draw Re-levels first
    }
    o = fopen("sqrt_map.eps","w");
    ado(o, 422.0, 422.0);
    std::fprintf(o,"211 211 translate\n 100 100 scale\n");
    draw_grid_eps(o, X, Y, M, N);
    // choose reasonable level list
    for (int k=-3;k<=3;k++) {
        std::fprintf(o,"0 0 0 RGB .01 W\n");
        // use a max_jump to block crossing the branch cut (approx)
        ContoArrayJump15(o, field, X, Y, M, N, (double)k, 2.0);
    }
    // Im levels
    for (int k=-3;k<=3;k++) {
        std::fprintf(o,"0 .6 0 RGB .01 W\n");
        ContoArrayJump15(o, field, X, Y, M, N, (double)k, 2.0);
    }
    std::fprintf(o,"showpage\n%%%%Trailer\n");
    fclose(o);

    // ---------- 3) log map (principal branch)
    for (int j=0;j<NY;j++) for (int i=0;i<NX;i++) {
        cd z(X[i], Y[j]);
        cd w = std::log(z); // principal log
        field[i + j*NX] = std::imag(w); // draw imaginary levels to show branch
    }
    o = fopen("log_map.eps","w");
    ado(o, 422.0, 422.0);
    std::fprintf(o,"211 211 translate\n 100 100 scale\n");
    draw_grid_eps(o, X, Y, M, N);
    // imaginary integer levels - show cut along negative real axis with jump suppression
    for (int k=-6;k<=6;k++) {
        if (k < 0) std::fprintf(o,"1 0 0 RGB .01 W\n"); else if (k>0) std::fprintf(o,"0 0 1 RGB .01 W\n"); else std::fprintf(o,".6 0 .6 RGB .02 W\n");
        ContoArrayJump15(o, field, X, Y, M, N, (double)k, 3.0); // skip edges with jump > 3
    }
    std::fprintf(o,"showpage\n%%%%Trailer\n");
    fclose(o);

    // cleanup
    free(X); free(Y); free(field);

    std::puts("Done: identity_map.eps, sqrt_map.eps, log_map.eps");
    return 0;
}
//</pre>