//
//+-------------------------------------------------------------+
//+ Mandelbrot Set via Continuous Potential Method with         +
//+ detection of Periods 1 & 2, January 30th 1994               +
//+ By Fausto Barbuto {BJ06@C53000.PETROBRAS.ANRJ.BR}, Brazil.  +
//+ Supports up to five screen types.                           +
//+ Needs SVGA256.BGI and SVGA256.H                             +
//+ REFERENCE:                                                  +
//+ Peitgen, H-O. & Saupe, D. - "The Science of Fractal Images" +
//+ Springer-Verlag, New York, 1988 ; Chap. 4, pp. 169-218      +
//+-------------------------------------------------------------+
//
#include <graphics.h>
#include <math.h>
#include <complex.h>
#include <stdio.h>
#include <conio.h>
#include "svga256.h"

int Vid; // Global variable

double MSetPot(double, double, int);

int huge DetectVGA256()
{
  printf("\nWhich video mode would you like to use? \n\n");
  printf(" 0 - 320x200x256\n");
  printf(" 1 - 640x400x256\n");
  printf(" 2 - 640x480x256\n");
  printf(" 3 - 800x600x256\n");
  printf(" 4 - 1024x768x256\n\n>");
  scanf("%d",&Vid);
  if((Vid<0) || (Vid)>4) Vid = 2;
  return Vid;
}

void main()
{
  int nx, ny, iy, ix, ipen, maxiter, iflag=0, iset;
  complex c;
  double xmin=-2.25, ymin=-1.25, xmax=0.75, ymax=1.25, cx, cy, potent;
  double diff=0.6482801, test1, test2;
  int graphdriver=DETECT, graphmode;

  clrscr();
  printf("\n Enter the number of allowed iterations (default=maximum=16000)\n\n>");
  scanf("%d",&maxiter);
  printf("\n Do you wish to plot Periods 1 & 2 with different colours?\n");
  printf("\n (1 = Yes ; <>1 = No) \n\n>");
  scanf("%d",&iset);
  if (iset != 1) iset = 0;
  if ((maxiter>=16000) || (maxiter<=0)) maxiter = 16000;

  installuserdriver("Svga256",DetectVGA256);
  initgraph(&graphdriver,&graphmode,"c:\\borlandc\\bgi");
  if (Vid == 0) { nx = 320; ny = 200; ymin = -1.250; ymax = 1.250;};
  if (Vid == 1) { nx = 640; ny = 400; ymin = -1.250; ymax = 1.250;};
  if (Vid == 2) { nx = 640; ny = 480; ymin = -1.125; ymax = 1.125;};
  if (Vid == 3) { nx = 800; ny = 600; ymin = -1.125; ymax = 1.125;};
  if (Vid == 4) { nx =1024; ny = 768; ymin = -1.125; ymax = 1.125;};
  for (iy=0;iy<=ny-1;iy++) {
    cy = ymin + iy*(ymax-ymin)/(ny-1);
    for (ix=0;ix<=nx-1;ix++) {
      cx = xmin + ix*(xmax-xmin)/(nx-1);
      c = complex(cx,cy);
//
//-------Checks limits for Period 1.
//
	 test1 = 2.0;
	 if ((cx >= -7.55e-1) && (cx <= 4.0e-1)) {
	   if ((cy >= -6.6e-1) && (cy <= 6.6e-1))
	      test1 = abs(1.0 - sqrt(1.0-4.0*c));
	 }
//
//-------Checks limits for Period 2.
//
	 test2 = 2.0;
	 if ((cx >= -1.275e0) && (cx <= -7.45e-1)) {
	   if ((cy >= -2.55e-1) && (cy <= 2.55e-1))
	      test2 = abs(4.0*(c+1.0));
	 }
//
      if (test1<=1.0) {
	potent = 0;
	iflag = 1;
	if (iset != 0) ipen = 126;
	else ipen = 32;
      }
      else if (test2<=1.0) {
	potent = 0;
	iflag = 1;
	if (iset != 0) ipen = 104;
	else ipen = 32;
      }
      else {
	potent = MSetPot(cx,cy,maxiter);
	iflag = 0;
      }
      if ((potent == 0.0) && (iflag==0))
	ipen = 32;
      else if ((potent !=0) && (iflag==0))
	ipen = (int)(33.0 + 15.0*(potent-33.0)/diff);
      putpixel(ix,iy,ipen);
    }
    if (kbhit()) break;
  }
  getch();
  closegraph();
}
double MSetPot(double cx, double cy, int maxiter)
{
   double x, y, x2, y2, temp, potential;
   int iter;

   x = cx;
   x2 = x*x;
   y = cy;
   y2 = y*y;
   iter = 0;

   do {
      temp = x2 - y2 + cx;
      y = 2.0*x*y + cy;
      x = temp;
      x2 = x*x;
      y2 = y*y;
      iter++;
   } while ((iter<maxiter) && ((x2+y2)<10000.0));
   if (iter<maxiter) potential = 0.5*log(x2+y2)/powl(2.0,iter);
   else potential = 0.0;
   return (potential);
}