分形曼德尔波集合图形

#include "glib.h"
#define KL 200
#define BOX 0.01
#define STEP 1
int sx,sy,dx=192,dy=192;
/*描绘区域*/
double x_min=-2.25,x_max=0.75,y_min=-1.5,y_max=1.5;
double x0,y0,a1,b1,a2,b2;

void coord(void)
{
  a1=(x_max-x_min)/(2*dx);
  b1=(x_max+x_min)/2;
  a2=(y_max-y_min)/(2*dy);
  b2=(y_max+y_min)/2;
}

void plot(int col)
{
  putpixel(sx,-sy,col);
  putpixel(sx,sy,col);
}

void mandelbrot(void)
{
  int k;
  double s;
  static double x[KL+1],y[KL+1];
  x[0]=x0;   y[0]=y0;
  for (k=1;k<=KL;k++)
  {
 x[k]=x[k-1]*x[k-1]-y[k-1]*y[k-1]+x0;  /*实部*/
 y[k]=2*x[k-1]*y[k-1]+y0;              /*虚部*/
 s=x[k]*x[k]+y[k]*y[k];
 if (s>=4.0)
   return;
  }
  for (k=1;k<=13;k++)    /*周期点数计算*/
  {
 if (fabs(x[KL]-x[KL-k])<BOX && fabs(y[KL]-y[KL-k])<BOX)
 {
   plot(k);
   return;
 }
  }
  plot(15);
}

void main()
{
  ginit();
  setpalette(8,14);
  set0(320,200);
  coord();
  /*坐标轴的描画*/
  setcolor(7);
  rectangle(-dx,-dy,dx,dy);
  for (sx=-dx+dx/6;sx<dx;sx+=dx/3)
  {
 setlinestyle(4,0x1111,1);
 line(sx,-dy,sx,dy);
  }
  for (sy=-dy+dy/3;sy<dy;sy+=dy/3)
  {
 setlinestyle(4,0x1111,1);
 line(-dx,sy,dx,sy);
  }
  for (sx=0;sx<=dx;sx+=STEP)
  {
 x0=a1*sx+b1;     /*x的初值*/
 for (sy=0;sy<=dy;sy+=STEP)
 {
   y0=a2*sy+b2;   /*y的初值*/
   mandelbrot();
 }
  }
  /*描画区域的显示*/
  printf("\nx_min =%51.f\n",x_min);
  printf("x_max =%5.1f\n",x_max);
  printf("y_min =%5.1f\n",y_min);
  printf("y_max =%5.1f\n",y_max);
  gend();
}