The Mandelbrot Set

The Mandelbrot Set

#include <stdio.h>
#include "dislin.h"

#define N 800
#define NITER 100

int iterate (double cx, double cy, int nmax);

float zmat[N][N];

main ()
{ float xscl[2], yscl[2], zscl[2];
  double cx, cy, xd, yd;
  int i, j, n, nmax = 100;

  metafl ("cons");
  scrmod ("revers");

  disini ();
  pagera ();
  hwfont ();

  axspos (300, 1900);
  ax3len (2200,1700,1700);

  xscl[0] = -2.f;
  xscl[1] = 1.f;
  yscl[0] = -1.f;
  yscl[1] = 1.f;
  zscl[0] = 1.f;
  zscl[1] = nmax;

  name ("Real parts", "x");
  name ("Imaginary parts", "y");
  name ("Number of iterations", "z");

  axsscl ("log", "z");
  labels ("log", "z");

  titlin ("The Mandelbrot Set", 4);

  setscl (xscl, 2, "x");
  setscl (yscl, 2, "y");
  setscl (zscl, 2, "z");

  graf3 (-2., 1., -2., 0.5, -1., 1., -1., 0.5, 0., 100., 10., 10.); 

  xd = (xscl[1] - xscl[0]) / (N - 1);
  yd = (yscl[1] - yscl[0]) / (N - 1);

  for (i = 0; i < N; i++)
  { cx = xscl[0] + i * xd;
    for (j = 0; j < N; j++)
    { cy = yscl[0] + j * yd;
      n = iterate (cx, cy, nmax);
      zmat[i][j] = n;
    }
  }

  crvmat ((float *) zmat, N, N, 1, 1);

  htitle (45);
  title ();
  disfin ();
}

int iterate (double cx, double cy, int nmax)
{ int n = 0;
  double x = 0., y = 0., x2 = 0., y2 = 0., zb = 0., zbmax = 4.;

  while (zb < zbmax && n < nmax) 
  { y = 2 * x * y + cy;
    x = x2 - y2 + cx;
    x2 = x * x;
    y2 = y * y;
    zb = x2 + y2;
    n++;
  }
  return n;
}

The Mandelbrot Set (Widget Program)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "dislin.h"

#define N 600

void myplot (int id);
void zoomplot (int id);
void resetplot (int id);
void undoplot (int id);
int  iterate (double cx, double cy);

float zmat[N][N], xscl[2], yscl[2], zscl[2], xold[2], yold[2], zold[N][N];
int iplot = 0, izoom = 0, iundo = 0,  niter = 100, izlog = 0, iold = 0;

int id_draw, id_xmin, id_xmax, id_ymin, id_ymax, id_zoom, id_reset, id_iter, 
    id_scale, id_progress, id_cancel, id_undo;

main ()
{ int ip, ip1, ip2, id;

  swgtit ("DISLIN Mandelbrot Plot");

  ip = wgini ("hori");
  swgwth (-70);
  ip1 = wgbas (ip, "vert");
  swgwth (-15);
  ip2 = wgbas (ip, "vert");

  swgdrw (2100./2970.);
  id_draw = wgdraw (ip1);

  id_xmin = wgltxt (ip2, "xmin:", "-2.000000", 50); 
  id_xmax = wgltxt (ip2, "xmax:", "1.000000", 50); 
  id_ymin = wgltxt (ip2, "ymin:", "-1.000000", 50); 
  id_ymax = wgltxt (ip2, "ymax:", "1.000000", 50); 

  id = wglab (ip2, " ");
  id_scale = wgbut (ip2, "Log. Colous", 0); 
  swgopt ("smooth", "pbar");

  id = wglab (ip2, " ");
  id = wglab (ip2, "Progress:");
  id_progress = wgpbar (ip2, 0., 100., 5.); 

  id = wglab (ip2, " ");
  id_iter = wgltxt (ip2, "Iterations:", "100", 40); 

  id = wglab (ip2, " ");
  id_zoom = wgpbut (ip2, "Zoom");
  swgcbk (id_zoom, zoomplot); 

  id_undo = wgpbut (ip2, "Undo Zoom");
  swgcbk (id_undo, undoplot); 

  id_cancel = wgpbut (ip2, "Cancel");
  id_reset = wgpbut (ip2, "Reset");
  swgcbk (id_reset, resetplot); 

  id = wglab (ip2, " ");
  id = wgquit (ip2);
  id = wgpbut (ip2, "Plot");
  swgcbk (id, myplot); 

  wgfin ();
}  

void myplot (int id)
{ int i, j, n, nclr, isel, nx1, ny1, nx2, ny2;
  double cx, cy, xd, yd;
  float xa = -1.0f, xe = 1.0f, xor = -1.0f, xstp = 0.2f, 
        ya = -1.0f, ye = 1.0f, yor = -1.0f, ystp = 0.2f,  
        za = 1.0f, ze = 100.0f, zor = 10.0f, zstp = 10.0f;  

  setxid (id_draw, "widget"); 

  xscl[0] = gwgflt (id_xmin);
  xscl[1] = gwgflt (id_xmax);
  yscl[0] = gwgflt (id_ymin);
  yscl[1] = gwgflt (id_ymax);
  niter = gwgint (id_iter);
  izlog = gwgbut (id_scale);  

  xd = (xscl[1] - xscl[0]) / (N - 1);
  yd = (yscl[1] - yscl[0]) / (N - 1);

  metafl ("cons");
  scrmod ("revers");

  disini ();
  if (izoom == 0) erase ();

  zscl[0] = 1.f;
  zscl[1] = (float) niter;

  setscl (xscl, 2, "x");
  setscl (yscl, 2, "y");
  setscl (zscl, 2, "z");

  nochek ();
  axspos (300, 1900);
  ax3len (2200,1700,1700);

  if (izlog == 1)
  { axsscl ("log", "z");
    labels ("log", "z");
  }
  else 
  { axsscl ("lin", "z");
    labels ("float", "z");
  }

  graf3 (xa, xe, xor, xstp, ya, ye, yor, ystp, za, ze, zor, zstp); 
  sendbf ();

  if (izoom == 0 && iundo == 0)  
  { for (i = 0; i < N; i++)
    { cx = xscl[0] + i * xd;
      doevnt ();
      if (gwgbut (id_cancel) == 1) 
      { swgbut (id_cancel, 0); 
        erase ();
        iplot = 0;
        disfin ();
        return;
      }        

      swgval (id_progress, i * 100.f / N);
      for (j = 0; j < N; j++)
      { cy = yscl[0] + j * yd;
        n = iterate (cx, cy);
        zmat[i][j] = n;
      }
    }
  }

  crvmat ((float *) zmat, N, N, 1, 1);

  if (izoom == 1)
  { csrrec (&nx1, &ny1, &nx2, &ny2);

    xold[0] = xscl[0];
    xold[1] = xscl[1];
    yold[0] = yscl[0];
    yold[1] = yscl[1];

    for (i = 0; i < N; i++)
      for (j = 0; j < N; j++)
	zold[i][j] = zmat[i][j];

    iold = 1;

    xscl[0] = xinvrs (nx1);
    xscl[1] = xinvrs (nx2);
    yscl[1] = yinvrs (ny1);
    yscl[0] = yinvrs (ny2);

    swgflt (id_xmin, xscl[0], 6);
    swgflt (id_xmax, xscl[1], 6);
    swgflt (id_ymin, yscl[0], 6);
    swgflt (id_ymax, yscl[1], 6);
  }

  disfin ();
  iplot = 1;
}

void zoomplot (int id)
{
  if (iplot == 0) return;      /* nothing to zoom */

  izoom = 1;
  myplot (id);                 /* just get the new scaling */
  if (iplot == 0) return;      /* maybe canceled */

  izoom = 0;
  myplot (id);                 /* replot with new scaling  */ 
}

void undoplot (int id)
{ int i, j;

  if (iold == 0) return;      /* nothing to do */

  xscl[0] = xold[0];
  xscl[1] = xold[1];
  yscl[0] = yold[0];
  yscl[1] = yold[1];

  swgflt (id_xmin, xscl[0], 6);
  swgflt (id_xmax, xscl[1], 6);
  swgflt (id_ymin, yscl[0], 6);
  swgflt (id_ymax, yscl[1], 6);

  for (i = 0; i < N; i++)
    for (j = 0; j < N; j++)
       zmat[i][j] = zold[i][j];

  iundo = 1;
  myplot (id);
  iundo = 0;
}

void resetplot (int id)
{
  xscl[0] = -2.f;
  xscl[1] = 1.f;
  yscl[0] = -1.f;
  yscl[1] = 1.f;

  swgflt (id_xmin, xscl[0], 6);
  swgflt (id_xmax, xscl[1], 6);
  swgflt (id_ymin, yscl[0], 6);
  swgflt (id_ymax, yscl[1], 6);

  myplot (id);
}

int iterate (double cx, double cy)
{ int n = 0;
  double x = 0., y = 0., x2, y2, zb = 0., zbmax = 4.;

  while (zb < zbmax && n < niter) 
  { x2 = x * x - y * y + cx;
    y2 = 2 * x * y + cy;
    x = x2;
    y = y2;
    n++;
    zb = x * x + y * y;
  }
  return n;
}

Rotation of a 3D Surface

#include <stdio.h>
#include <math.h>
#include "dislin.h"

#define N 50

static int id_draw, id_but, id_lab1, id_lab2, id_scl1, id_scl2;
static float xrot1 = 50.0f, xrot2 = 40.0f;
static float xray[N], yray[N], zmat[N][N];

void myplot (int id);

int main()
{ int ip, ip1, ip2, i, j;
  double fpi=3.1415927/180., step, x, y;

  step = 360./(N-1);
  for (i = 0; i < N; i++)
  { x = i*step;
    xray[i] = x;
    for (j = 0; j < N; j++)
    { y = j*step;
      yray[j] = y;
      zmat[i][j] = (float) (2*sin(x*fpi)*sin(y*fpi));
    }
  }

  swgtit ("DISLIN 3D Plot");
  ip = wgini ("hori");
  swgwth (-50);
  ip1 = wgbas (ip, "vert");
  swgwth (-15);
  ip2 = wgbas (ip, "vert");

  swgdrw ((float) (2100./2970.));
  id_lab2 = wglab (ip1, "DISLIN 3D Plot:");
  id_draw = wgdraw (ip1);

  swgopt ("track", "scroll");
  id_lab1 = wglab (ip2, "Azimuthal Rotation Angle:");
  id_scl1 = wgscl (ip2, " ", 0.f, 360.f, 50.f, -1);
  id_lab1 = wglab (ip2, "Polar Rotation Angle:");
  id_scl2 = wgscl (ip2, " ", -90.f, 90.f, 40.f, -1);
  swgcbk (id_scl1, myplot); 
  swgcbk (id_scl2, myplot); 

  myplot (id_scl1); 
  wgfin ();
  return 0;
}  

void myplot (int id)
{
  if (id == id_scl1)
    xrot1 = gwgscl (id);
  else if (id == id_scl2)
    xrot2 = gwgscl (id);

  scrmod ("revers");
  setxid (id_draw, "widget");
  metafl ("xwin");
  disini ();
  erase();
  hwfont ();

  name ("X-axis", "x");
  name ("Y-axis", "y");
  name ("Z-axis", "z");
  labdig (-1, "xyz");

  axspos (585, 1800);
  axslen (1800, 1800);

  view3d (xrot1, xrot2, 6.0f, "angle");
  height (40);
  graf3d(0.f,360.f,0.f,90.f,0.f,360.f,0.f,90.f,-3.f,3.f,-3.f,1.f);
  surshd(xray, N, yray, N, (float *) zmat);
  disfin ();
}