//
//+-------------------------------------------------------------------+
//+ PROGRAMA PMAPS.CPP                                                +
//+ By Fausto A. A. Barbuto - Rio de Janeiro, BRAZIL, March 22, 1994. +
//+ E-mail: BJ06@C53000.PETROBRAS.ANRJ.BR                             +
//+                                                                   +
//+ Plots rational maps from complex functions of the form            +
//+ h(z,c) = [f(z,c)/g(z,c)]^2. 12 functions are proposed.            +
//+ Uses SVGA256.BGI and supports up to five video screens.           +
//+ I henceforward challenge the future users to create their own     +
//+ rational funcions.                                                +
//+-------------------------------------------------------------------+
//
#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 *);

int Vid;  //Global variable

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()
{
     double pmin, pmax, qmin, qmax;
     double ypy,x,y,xp,yp,p,q,r;
     double deltap, deltaq;
     register int npix, npiy, k, np, nq, npy, ipen;
     int index, istart;
     unsigned int maxiter;
     complex c, z, za, z2;
     int graphdriver=DETECT, graphmode;

     clrscr();
     printf("\n                Select a formulation:   \n\n");
     printf(" 1: ((z^2+c-1)/(2z+c-2))^2      2: ((z^2+c+1)/(z+c-2))^2\n");
     printf(" 3: ((z^2+c)/(z-c))^2           4: ((z+c)/(z-c))^2\n");
     printf(" 5: ((z^3-2z-c)/(2z+c))^2       6: ((z^2+c-2)/(2z+c+1))^2\n");
     printf(" 7: ((z^2+c-2)/(z+c+1))^2       8: ((z^2+2c+5)/(z^2+c))^2\n");
     printf(" 9: ((z^3+2c+9)/(z^2-3z-c))^2  10: ((z^2+2c-3)/(2z-c))^2\n");
     printf("11: ((z^2+2(c-1))/(z-c))^2     12: ((z^3+z^2-c)/(z-c))^2\n");
     printf("\n\n>");
     scanf("%d",&index);
     printf("\n  Number of iterations (maximum=default=16384)?");
     printf("\n  (suggested : 64 or higher)\n\n>");
     scanf("%d",&maxiter);
     printf("\n  Set starting colour (>=0 , <=256)?");
     printf("\n  (suggested: 33, 16, 40, 145)\n\n>");
     scanf("%d",&istart);
     if((maxiter>16384) || (maxiter)<=0) maxiter = 16384;
     if((istart)<=0 || (istart>256)) istart = 0;
     clrscr();

     installuserdriver("Svga256",DetectVGA256);
     initgraph(&graphdriver, &graphmode, "c:\\borlandc\\bgi");

     if (Vid == 0) { npix = 320; npiy = 200;}
     if (Vid == 1) { npix = 640; npiy = 400;}
     if (Vid == 2) { npix = 640; npiy = 480;}
     if (Vid == 3) { npix = 800; npiy = 600;}
     if (Vid == 4) { npix =1024; npiy = 768;}
     if ((Vid>4) || (Vid<0)) { npix = 640; npiy = 480;}

//---Sets plotting window for each rational function.

     if(index == 1) {pmin=-0.50; pmax=4.00;  qmin=-2.25; qmax=2.25;}
     if(index == 2) {pmin=-0.50; pmax=5.50;  qmin=-3.00; qmax=3.00;}
     if(index == 3) {pmin=-0.50; pmax=5.50;  qmin=-3.00; qmax=3.00;}
     if(index == 4) {pmin=-0.50; pmax=5.50;  qmin=-3.00; qmax=3.00;}
     if(index == 5) {pmin=-6.00; pmax=0.00;  qmin=-3.00; qmax=3.00;}
     if(index == 6) {pmin=-5.60; pmax=-5.5;  qmin=-0.04; qmax=0.04;}
     if(index == 7) {pmin=-8.25; pmax=-5.25; qmin=-1.50; qmax=1.50;}
     if(index == 8) {pmin=-2.00; pmax=-1.00; qmin=-0.50; qmax=0.50;}
     if(index == 9) {pmin=-14.0; pmax=-6.00; qmin=-4.00; qmax=4.00;}
     if(index ==10) {pmin=-6.00; pmax=-4.00; qmin=-0.80; qmax=0.80;}
     if(index ==11) {pmin=-6.00; pmax=-5.00; qmin=-0.50; qmax=0.50;}
     if(index ==12) {pmin=-12.0; pmax=-8.00; qmin=-2.00; qmax=2.00;}

     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);
//
	 do {
	   z = complex(x,y);

//         ***** Definition of the Functions *****

//---------See "The Science of Fractal Images", Fig. 4.26a, pg. 213
//---------I am sorry, but this is difficult to match at the screen
//---------without distortion. I suggest to use maxiter = 50.
	   if(index== 1) za = (z*z + c - 1.0)/(2*(z - 1.0) + c);
//
	   if(index== 2) za = (z*z + c + 1.0)/(z + c - 2.0);
	   if(index== 3) za = (z*z + c)/(z - c);
//
//---------Dust.
	   if(index== 4) za = (z + c)/(z - c);
//
//---------Arabesque.
	   if(index== 5) za = (z*z*z - 2.0*z - c)/(2.0*z + c);
//
//---------What the hell is the Mandelbrot Set doing here???
	   if(index== 6) za = (z*z + c - 2.0)/(2.0*z + c + 1.0);
//
	   if(index== 7) za = (z*z + c - 2.0)/(z + c + 1.0);
	   if(index== 8) {z2=z*z; za = (z2 + 2.0*c + 5.0)/(z2 + c);}
	   if(index== 9) {z2=z*z; za = (z*z2 + 2.*c + 9.)/(z2 - 3.*z - c);}
//
//---------A Mandelbrot Set? Again??
	   if(index==10) za = (z*z + 2.0*c - 3.0)/(2.0*z - c);
//
//---------A variation on Formula 10 shows another... guess what?
	   if(index==11) za = (z*z + 2.0*(c - 1.0))/(z - c);
//
	   if(index==12) za = (z*z*z + z*z - c)/(z - c);

	   z = za*za;
	   r = abs(z);
	   k++;
//
//---------Plots escaping points.
//
	   if (r >= maxiter) {
	     ipen = istart + k;
	     if (qmin == -qmax) {
	       putpixel(np,nq,ipen);
	       putpixel(np,npiy-nq-1.0,ipen);
	     }
	     else
	       putpixel(np,nq,ipen);
	   }
//
//---------Plots atracted points.
//
	   if (k == maxiter) {
	     ipen = 1;
	    if (qmin == -qmax) {
	      ypy = double(npiy) - nq - 0.5;
	      putpixel(np,ypy,ipen);
	      putpixel(np,nq,ipen);
	    }
	    else
	      putpixel(np,nq,ipen);
	   }
//
//---------Prepares a new iteration.
//
	   x = real(z);
	   y = imag(z);
	 } while (r <= maxiter && k<=maxiter);
       }
       if(kbhit()) break;
     }
//
//---Clean up and end.
//
     getch();
     closegraph();
}