Difference between revisions of "File:LogTestMapChatGPTv12.jpg"

This is 12th attempt to use
version 5 of conto02.cin
to generate the contour plot of function ln
handling both the cutline and the singularity.

The mesh points are intentionally misplaced from their native position in order to see defects of the implementation.

The picture is loaded for use in article LogTestMaps

05.cin

// conto02.cin  -- portable array-based contour tracer (C++98 compatible)
//
// Usage:
//   // compute arrays: X[0..M], Y[0..N], field[ (M+1)*(N+1) ] indexed as i + j*(M+1)
//   // set color/linewidth in caller's EPS stream (e.g., "/W { setlinewidth } def" etc.)
//   conto_array(FILE* out, double *field, double *X, double *Y, int M, int N, double level);
//

#ifndef CONTO02_CIN_ONCE
#define CONTO02_CIN_ONCE

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

static const double EPS_ZERO = 1e-14;

struct Pt { double x, y; };

// safe interpolation between two edge endpoints (A->B order)
static inline Pt interp_edge(double xA,double yA,double fA,
                             double xB,double yB,double fB,double level)
{
    Pt p;
    double denom = (fB - fA);
    if (fabs(denom) < EPS_ZERO) { p.x = 0.5*(xA + xB); p.y = 0.5*(yA + yB); }
    else {
        double t = (level - fA) / denom;
        p.x = xA + t*(xB - xA);
        p.y = yA + t*(yB - yA);
    }
    return p;
}

// bilinear discriminant (unused now, kept for future option)
static inline double amb_decider(double f00,double f10,double f01,double f11){
    return f00 * f11 - f10 * f01;
}

// Draw segment to EPS (caller must ensure PostScript macros exist)
static inline void draw_seg(FILE* out, const Pt &a, const Pt &b) {
    std::fprintf(out, "newpath %.12g %.12g moveto %.12g %.12g lineto stroke\n", a.x, a.y, b.x, b.y);
}

// Main array-based function
void conto_array(FILE* out,
                 double *field, double *X, double *Y,
                 int M, int N, double level)
{
    if(!out || !field || !X || !Y) return;
    if(M <= 0 || N <= 0) return;

    int NX = M + 1;
    int NY = N + 1;

    for (int j = 0; j < N; ++j) {
        double y0 = Y[j], y1 = Y[j+1];
        for (int i = 0; i < M; ++i) {
            double x0 = X[i], x1 = X[i+1];

            // 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];

            // fast skip (all strictly same side)
            if ((f00 - level) > 0.0 && (f10 - level) > 0.0 && (f11 - level) > 0.0 && (f01 - level) > 0.0) continue;
            if ((f00 - level) < 0.0 && (f10 - level) < 0.0 && (f11 - level) < 0.0 && (f01 - level) < 0.0) continue;

            // compute edge intersections with consistent endpoint ordering:
            // edge 0 = bottom (x0,y0)->(x1,y0)
            // edge 1 = right  (x1,y0)->(x1,y1)
            // edge 2 = top    (x0,y1)->(x1,y1)
            // edge 3 = left   (x0,y0)->(x0,y1)
            Pt edge[4];
            bool has[4] = {false,false,false,false};

            if ((f00 - level) * (f10 - level) < 0.0) {
                edge[0] = interp_edge(x0,y0,f00, x1,y0,f10, level); has[0] = true;
            }
            if ((f10 - level) * (f11 - level) < 0.0) {
                edge[1] = interp_edge(x1,y0,f10, x1,y1,f11, level); has[1] = true;
            }
            if ((f11 - level) * (f01 - level) < 0.0) {
                edge[2] = interp_edge(x0,y1,f01, x1,y1,f11, level); has[2] = true;
            }
            if ((f01 - level) * (f00 - level) < 0.0) {
                edge[3] = interp_edge(x0,y0,f00, x0,y1,f01, level); has[3] = true;
            }

            // collect present edges in order bottom,right,top,left
            int idxs[4]; int k = 0;
            for (int e = 0; e < 4; ++e) if (has[e]) idxs[k++] = e;

            if (k == 2) {
                // straightforward case
                Pt a = edge[idxs[0]], b = edge[idxs[1]];
                draw_seg(out, a, b);
                continue;
            }

            if (k == 4) {
                // Ambiguous saddle: default behavior -> draw both diagonals (makes an X crossing)
                // Diagonals: (bottom edge midpoint -> top edge midpoint)
                //            (left edge midpoint   -> right edge midpoint)
                Pt pb = edge[0], pr = edge[1], pt = edge[2], pl = edge[3];
                draw_seg(out, pb, pt); // bottom -> top
                draw_seg(out, pl, pr); // left -> right
                continue;
            }

            // handle equality corner cases (within EPS_ZERO)
            bool c00eq = fabs(f00 - level) < EPS_ZERO;
            bool c10eq = fabs(f10 - level) < EPS_ZERO;
            bool c01eq = fabs(f01 - level) < EPS_ZERO;
            bool c11eq = fabs(f11 - level) < EPS_ZERO;

            // opposite corners equal -> draw full diagonal
            if ((c00eq && c11eq) || (c10eq && c01eq)) {
                Pt a,b;
                if (c00eq && c11eq) { a.x = x0; a.y = y0; b.x = x1; b.y = y1; }
                else { a.x = x1; a.y = y0; b.x = x0; b.y = y1; }
                draw_seg(out, a, b);
                continue;
            }

            // single corner equals -> do bishop move small diagonal
            if ( (c00eq && !c10eq && !c01eq && !c11eq) ||
                 (c10eq && !c00eq && !c11eq && !c01eq) ||
                 (c11eq && !c10eq && !c01eq && !c00eq) ||
                 (c01eq && !c00eq && !c11eq && !c10eq) ) {
                Pt p,q;
                if (c00eq) { p.x = x0; p.y = y0; q.x = x1; q.y = y1; }
                else if (c10eq) { p.x = x1; p.y = y0; q.x = x0; q.y = y1; }
                else if (c11eq) { p.x = x1; p.y = y1; q.x = x0; q.y = y0; }
                else { p.x = x0; p.y = y1; q.x = x1; q.y = y0; }
                draw_seg(out, p, q);
                continue;
            }

            // Otherwise: degenerate/no-draw
        }
    }
}

Generator

#include <stdio.h>
#include <stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
#include<complex>
typedef std::complex<double> z_type;
#define Re(x) x.real()
#define Im(x) x.imag()
#define I z_type(0.,1.)
#include "ado.cin"
#include "05.cin"

z_type F(z_type z) {return log(z);}

int main(){ int j,k,m,n; DB x,y, p,q, t; z_type z,c,d, cu,cd;
int M=41,M1=M+1;
int N=41,N1=N+1;
DB X[M1],Y[N1]; DB *g, *f, *w; // w is working array.
g=(DB *)malloc((size_t)((M1*N1)*sizeof(DB)));
f=(DB *)malloc((size_t)((M1*N1)*sizeof(DB)));
w=(DB *)malloc((size_t)((M1*N1)*sizeof(DB)));
char v[M1*N1]; // v is working array
//FILE *o;o=fopen("QuadraticMapByChatGPTtest11.eps","w");  ado(o,422,422);
//FILE *o;o=fopen("F1zChatGPTtest.eps","w");  ado(o,422,422);
//FILE *o;o=fopen("09.eps","w");  ado(o,422,422);
FILE *o;o=fopen("12.eps","w");  ado(o,422,422);
fprintf(o,"210 210 translate\n 100 100 scale\n");
fprintf(o,"1 setlinejoin 2 setlinecap\n");
fprintf(o,"/W { setlinewidth } def\n");
fprintf(o,"/RGB { setrgbcolor } def\n");
fprintf(o,"/S { stroke } def\n");
#define M(x,y) fprintf(o,"%6.3f %6.3f M\n",0.+x,0.+y); 
#define L(x,y) fprintf(o,"%6.3f %6.3f L\n",0.+x,0.+y); 
DO(m,M1) X[m]=-2+.1*(m-.4);
DO(n,N1) Y[n]=-2+.1*(n-.3); 
for(m=-2;m<3;m++){M(m,-2) L(m,2)}
for(n=-2;n<3;n++){M(-2,n) L(2,n)}
 fprintf(o,".01 W 0 0 0 RGB S\n");
DO(m,M1)DO(n,N1){      g[m*N1+n]=999;
                       f[m*N1+n]=999;}
DO(m,M1){x=X[m]; if(m/10*10==m) printf("x=%6.3f\n",x);
DO(n,N1){y=Y[n]; z=z_type(x,y); //if(abs(z+2.)>.019)
  c=F(z);
 p=Re(c); q=Im(c); if(p>-99 && p<99 )
//{ g[m*N1+n]=p;f[m*N1+n]=q;}
{ g[m+N1*n]=p;f[m+M1*n]=q;}
       }}


// handdle the cut line and the singularity
DO(m,M)DO(n,N1){ if( abs(f[m+N1*n]-f[m+1+N1*n])>2. || abs(g[m+N1*n]-g[m+1+N1*n])>1) {f[m+N1*n]=f[m+1+N1*n]=g[m+N1*n]=g[m+1+N1*n]=NAN;}  }
DO(m,M1)DO(n,N){ if( abs(f[m+N1*n]-f[m+N1*(n+1)])>2. || abs(g[m+N1*n]-g[m+N1*(n+1)])>1) {f[m+N1*n]=f[m+1+N1*n]=g[m+N1*n]=g[m+N1*(n+1)]=NAN;}}

fprintf(o,"1 setlinejoin 1 setlinecap\n"); p=16.;q=8.;
fprintf(o,".01 W 0 .6 0 RGB\n"); for(m=-3;m<3;m++)for(n=2;n<10;n+=2)conto_array(o,f,X,Y,M,N,(m+.1*n));
fprintf(o,".01 W .8 0 0 RGB\n"); for(m=0;m<3;m++) for(n=2;n<10;n+=2)conto_array(o,g,X,Y,M,N,-(m+.1*n));
fprintf(o,".01 W 0 0 .8 RGB\n"); for(m=0;m<3;m++) for(n=2;n<10;n+=2)conto_array(o,g,X,Y,M,N, (m+.1*n));
fprintf(o,".02 W 0 0 .8 RGB\n"); for(m= 1;m<5;m++) conto_array(o,f,X,Y,M,N, (0.+m)); 
fprintf(o,".02 W .8 0 0 RGB\n"); for(m= 1;m<5;m++) conto_array(o,f,X,Y,M,N, (0.-m));
fprintf(o,".02 W .5 0 .5 RGB\n");                   conto_array(o,f,X,Y,M,N, (0.  ));
fprintf(o,".02 W 0 0 0 RGB\n"); for(m=-10;m<11;m++)conto_array(o,g,X,Y,M,N,(0.+m));
fprintf(o,"showpage\n"); fprintf(o,"%c%cTrailer\n",'%','%');
fclose(o);  free(f); free(g); free(w);
      system("epstopdf 12.eps"); 
      system(    "open 12.pdf"); //for macintosh
}

References

Keywords

«Ado.cin», «ChatGPT», «Conto.cin», «Contour plot», «Log», «LogTestMaps‎‎», «ContoTestMaps»,