#include <stdio.h>
#include <conio.h>
#include <graphics.h>
#include <math.h>
#include <complex.h>
#include "svga256.h"
//void far initgraph(int far *,int far *,char far *);
//    PROGRAMA MANDELBROT_SHADINGS_1
//
//    By Fausto A. A. Barbuto, December 1993
//    (After a program by Fausto Barbuto and Jay R. Hill, 1993)
//
//    Plots the Mandelbrot Set with period checking.
//
//    Period 1 (also known as Main Cardioid) ---> is in colour
//    shadings from green to blue;
//
//    Formulation:
//               abs(1.0 - sqrt(1.0-4.0*c)) <= 1.0
//
//    Period 2 (also known as Main Circle) ---> is in colour
//    shadings from blue to orange;
//
//    Formulation:
//               abs(4.0*(c+1.0)) <= 1.0
//
//    Period 3 Buds (no generic name: periodic detected by Fausto Barbuto
//                   and Jay R. Hill, 1993) ---> is in dark blue;
//    Formulation:
//               p1*p1 + q1*q1 <= 0.0089138
//               p1 = p + 0.12486
//               q1 = fabs(q) - 0.7439
//
//    Period 4 (no generic name: periodic detected by Fausto Barbuto, 1993)
//             ---> is in dark blue;
//    Formulation:
//               abs(c+1.309) <= 0.0588
//
//    Other MSet points which do not belong to any of these periods are
//    coloured in dark blue.
//
//    Escaped points (external to the MSet) are in shadings from blue to
//    yellow.
//
//    PS: Needs SVGA256.BGI Graphics Driver and SVGA screen with 1024 x 768
//        points and 256 colours.
//
 
int huge DetectVGA256()
{
  int Vid=4;
  return Vid;
}
 
void main()
{
      double pmin=-2.05, pmax=0.5, qmin=-1.125, qmax=1.125, fact=1.0;
      double ypy,x,y,xp,yp,p,q,p1,q1,ya,xkp1,ykp1,r,c_scr=0.833333;
      double deltap, deltaq, r1=8.9138e-3, r2=5.88e-2;
      double test1, test2, test3, test4;
      register int npix=1024, npiy=768;
      register long int maxiter=512;
      register int k, np, nq, npy, ipen, icolour;
      complex c;
      int graphdriver=DETECT, graphmode;
 
      installuserdriver("Svga256",DetectVGA256);
      initgraph(&graphdriver, &graphmode, "c:\\borlandc\\bgi");
      if(fact>=1.0 || fact <=0.0)
        fact = 1.0;
      else {
        npix = (int)(npix*fact);
        npiy = (int)(npiy*fact);
      }
      ypy = (double)npiy - 0.5;
      deltap = (pmax-pmin)/(npix-1);
      deltaq = (qmax-qmin)/(npiy-1);
 
      if(qmin==-qmax)
         npy = npiy/2;
      else
         npy = npiy;
 
     cleardevice();
     for (np=0; np<=npix-1; np++) {
       p = pmin + (double)np*deltap;
       for (nq=0; nq<=npy-1; nq++) {
         q = qmin + (double)nq*deltaq;
         k  = 0;
         x = 0.0;
         y = 0.0;
         c = complex(p,q);
//
//-------Checks limits for Period 1.
//
         test1 = 2.0;
         if ((p >= -7.55e-1) && (p <= 4.0e-1)) {
           if ((q >= -6.6e-1) && (q <= 6.6e-1))
              test1 = abs(1.0 - sqrt(1.0-4.0*c));
         }
//
//-------Checks limits for Period 2.
//
         test2 = 2.0;
         if ((p >= -1.255e0) && (p <= -7.45e-1)) {
           if ((q >= -2.55e-1) && (q <= 2.55e-1)) test2 = abs(4.0*(c+1.0));
         }
//
//-------Checks limits for Period Bud 3.
//
         test3 = 2.0;
         if ((p >= -2.4e-1) && (p <= 0.0e0)) {
           if ((q >= -8.5e-1) && (q <= 8.5e-1)) {
             p1 = p + 1.2486e-1;
             q1 = fabs(q) - 7.4396e-1;
             test3 = p1*p1 + q1*q1;
           }
         }
//
//-------Checks limits for Period 4.
//
         test4 = 2.0;
         if ((p >= -1.37e0) && (p <= -1.245e0)) {
           if ((q >= -6.1e-2) && (q <= 6.1e-2)) test4 = abs(c+1.309);
         }
//
         if((test1<=1.0) || (test2<=1.0) || (test3<=r1) || (test4<=r2)) {
            if (test3<=r1)  ipen=1;
            if (test4<=r2)  ipen=1;
            xp = c_scr*(double)np;
            yp = (double)nq;
            if (test1<=1.0) {
              ipen = 33 + (int)(15.0*test1);
            }
            if (test2<=1.0) {
              ipen = 42 - (int)(10.0*test2);
            }
            if (qmin == -qmax) {
              putpixel(xp,yp,ipen);
              putpixel(xp,npiy-yp-1.0,ipen);
            }
            else
              putpixel(xp,yp,ipen);
            }
         else {
            do {
             xkp1 = (x+y)*(x-y) + p;
             ya   = x*y;
             ykp1 = ya + ya + q;
             r    = xkp1*xkp1 + ykp1*ykp1;
             k++;
//
//    If R > M, points escape towards infinity.
//    Se R > M, os pontos escapam para o infinito.
//
             if (r >= maxiter) {
               ipen = 28 + k;
               xp = c_scr*(double)np;
               yp = (double)nq;
               if (qmin == -qmax) {
                 putpixel(xp,yp,ipen);
                 putpixel(xp,npiy-yp-1.0,ipen);
               }
               else
                 putpixel(xp,yp,ipen);
             }
//
//    Points which do not belong to any period: Colour=1=Blue
//    Pontos que nao pertencem a quaisquer periodos: coloracao negra.
//
             if (k == maxiter) {
               xp = c_scr*(double)np;
               yp = (double)nq;
               ipen = 1;
               if (qmin == -qmax) {
                 ypy = double(npiy) - yp - 0.5;
                 putpixel(xp,ypy,ipen);
                 putpixel(xp,yp,ipen);
               }
               else
                 putpixel(xp,yp,ipen);
             }
//
//    Returns if no convergence is achieved on either escape or atraction.
//    Retorna se nao houver convergencia na fuga ou atracao.
//
             x = xkp1;
             y = ykp1;
           } while (r <= maxiter && k<=maxiter);
         }
       }
       if(kbhit()) break;
     }
     getch();
     closegraph();
}