//
//+-----------------------------------------------------------------+
//+ Program MAGNETS.CPP                                             +
//+ Plots 18 magnet-like fractals in the x-plane for different q's. +
//+ (16 for Model I, 2 for Model II)                                +
//+                                                                 +
//+ By: Fausto Arinos de Almeida Barbuto                            +
//+ E-mail: BJ06@C53000.PETROBRAS.ANRJ.BR                           +
//+ Rio de Janeiro, BRAZIL, April 4th 1994                          +
//+                                                                 +
//+ Needs SVGA256.H and SVGA256.BGI, supports up to five screens.   +
//+                                                                 +
//+ REFERENCE: Peitgen, H.-O. and Richter, P.H.:                    +
//+            "The Beauty of Fractals", Springer-Verlag, 1986,     +
//+            pg. 194.                                             +
//+-----------------------------------------------------------------+
//
#include <stdio.h>
#include "svga256.h"
#include <conio.h>
#include <graphics.h>
#include <math.h>
#include <complex.h>

//void far initgraph(int far *,int far *,char far *);

int Vid; // Global variable //

int huge DetectVGA256()
{
  printf("\n Which 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 xmin, xmax, ymin, ymax, fact=1.0;
      double ypy, x, y, r, deltap, deltaq;
      int istart, maxiter, Model;
      register int npix, npiy, k, np, nq, npy, ipen;
      complex c, z, q, a, b, ab, z2;
      int graphdriver=DETECT, graphmode, index;

      clrscr();
      printf("\n  MAGNET-LIKE MODELS");
      printf("\n  Reference: Peitgen, H.-O., and Richter, P.H. :");
      printf("\n  _The Beauty of Fractals_, Springer-Verlag, 1986, pg. 194");
      printf("\n\n  By Fausto A. A. Barbuto, April 4, 1994");
      printf("\n  E-mail: BJ06@C53000.PETROBRAS.ANRJ.BR \n\n\n");
      printf("\n  SELECT A MODEL TYPE: \n");
      printf("\n  Model I:");
      printf("\n     z-> ((z^2 + q - 1)/(2z + q - 2))^2 \n");
      printf("\n  Model II:");
      printf("\n     z-> ((z^3 + 3az + ab)/(3z^2 + 3bz + q^2 - 3q + 3))^2");
      printf("\n     (where a = q-1 and b = q-2) \n\n");
      printf("\n  (Model I = 1 ; Model II = 2)\n\n>");
      scanf("%d",&Model);
      clrscr();
      if ((Model != 1) && (Model != 2)) Model = 1;
      if (Model == 1) {
       printf("\n             Model I: Select a formulation:   \n\n\n");
       printf("   1: q = -1.0       -6<=Re(z)<=4       -5<=Im(z)<=5 \n");
       printf("   2: q = -0.1       -6<=Re(z)<=4       -5<=Im(z)<=5 \n");
       printf("   3: q =  0.0       -6<=Re(z)<=4       -5<=Im(z)<=5 \n");
       printf("   4: q =  1.0       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("   5: q =  1.2       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("   6: q =  1.6       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("   7: q =  2.0       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("   8: q =  2.5       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("   9: q =  2.9       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("  10: q =  3.0       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("  11: q =  3.1       -5<=Re(z)<=5       -5<=Im(z)<=5 \n");
       printf("  12: q =  4.0       -4<=Re(z)<=6       -5<=Im(z)<=5 \n");
       printf("  13: q =  4.0     -5.5<=Re(z)<=7.9     -5<=Im(z)<=5\n");
       printf("  14: q = -0.1      0.8<=Re(z)<=2.0  2.635<=Im(z)<=3.535\n");
       printf("  15: q=~1.1+2.07i -0.5<=Re(z)<=1.5  -1.45<=Im(z)<=0.7\n");
       printf("  16: q=1.21+0.01i -2.1<=Re(z)<=-0.3  -0.5<=Im(z)<=0.5625");
       printf("\n\n>");
       scanf("%d",&index);
       clrscr();
      }
      else
      if (Model == 2) {
       printf("\n             Model II: Select a formulation:   \n\n\n");
       printf("   1: q = 2.0       -15<=Re(z)<=15     -11<=Im(z)<=11\n");
       printf("   2: q = 1.2+2i   -3.2<=Re(z)<=3.8   -3.3<=Im(z)<=1.9\n>");
       scanf("%d",&index);
      }
      if ((index>16) || (index<1)) index = 1;
      clrscr();
      printf("\n Number of iterations (maximum=16000, default=256)?");
      printf("\n (Suggested: 64 or higher) \n\n>");
      scanf("%d",&maxiter);
      printf("\n Set starting colour (>=0 , <=256)?");
      printf("\n (Suggested: 33, 30, 16, 40, 145)\n\n>");
      scanf("%d",&istart);
      if((maxiter>16000) || (maxiter)<=0) maxiter = 256;
      if((istart)<=0 || (istart>256)) istart = 33;

      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;}
//
//    Defaults for the plotting window.
//
      xmin = -5.0; xmax = 5.0; ymin = -5.0; ymax = 5.0;
//
      if (Model == 1) {
	if (index == 1) {xmin=-6.0, xmax=4.0; q=complex(-1.0,0.0);}
	else
	if (index == 2) {xmin=-6.0, xmax=4.0; q=complex(-0.1,0.0);}
	else
	if (index == 3) {xmin=-6.0, xmax=4.0; q=complex(0.0,0.0);}
	else
	if (index == 4) q=complex(1.0,0.0);
	else
	if (index == 5) q=complex(1.2,0.0);
	else
	if (index == 6) q=complex(1.6,0.0);
	else
	if (index == 7) q=complex(2.0,0.0);
	else
	if (index == 8) q=complex(2.5,0.0);
	else
	if (index == 9) q=complex(2.9,0.0);
	else
	if (index ==10) q=complex(3.0,0.0);
	else
	if (index ==11) q=complex(3.1,0.0);
	else
	if (index ==12) {xmin=-4.0, xmax=6.0; q=complex(4.0,0.0);}
	else
	if (index ==13) {xmin=-5.5, xmax=7.9; q=complex(4.0,0.0);}
	else
	if (index ==14) {xmin=0.8, xmax=2.0; ymin=2.635; ymax=3.535;
			 q=complex(-0.1,0.0);}
	else
	if (index ==15) {xmin=-0.5, xmax=1.5; ymin=-1.45; ymax=0.7;
			 q=complex(1.09582,2.07142);}
	else
	if (index ==16) {xmin=-2.1, xmax=-0.3; ymin=-0.5; ymax=0.5625;
			 q=complex(1.21,0.01);}
	else {q=complex(1.2,0.0);}
      }
      else if (Model == 2) {
	if (index ==1) {xmin=-15., xmax=15.; ymin=-11.0; ymax=11.0;
			q=complex(2.0,0.0);}
	else
	if (index ==2) {xmin=-3.2, xmax=3.8; ymin=-3.3; ymax=1.9;
			q=complex(1.20,2.0);}
	else {xmin=-3.2, xmax=3.8; ymin=-3.3; ymax=1.9;
	      q=complex(1.20,2.0);}
      }

      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 = (xmax-xmin)/(npix-1);
      deltaq = (ymax-ymin)/(npiy-1);

      if(ymin==-ymax)
	 npy = npiy/2;
      else
	 npy = npiy;

     cleardevice();
     for (np=0; np<=npix-1; np++) {
       x = xmin + (double)np*deltap;
       for (nq=0; nq<=npy-1; nq++) {
	 y = ymin + (double)nq*deltaq;
	 k  = 0;
	 z = complex(x,y);
	 do {
//
//         Definition for Model I magnets.
//
	   if (Model == 1) z = (z*z + q - 1.0)/(2.0*z + q - 2.0);
	   else
//
//         Definition for Model II magnets.
//
	   if (Model == 2) {
	     a = q - 1.0;
	     b = q - 2.0;
	     ab = a*b;
	     z2 = z*z;
	     z = (z*z2 + 3.0*a*z + ab)/(3.0*z2 + 3.0*b*z + ab + 1.0);
	   }
	   z = z*z;
	   r = abs(z);
	   k++;
//
//         Points which escape towards infinity.
//
	   if (r >= maxiter) {
	     ipen = istart + k;
	     if (ymin == -ymax) {
	       putpixel(np,nq,ipen);
	       putpixel(np,npiy-nq-1,ipen);
	     }
	     else
	       putpixel(np,nq,ipen);
	   }
//
//         Converging points.
//
	   if (k == maxiter) {
	     ipen = 33;
	     if (ymin == -ymax) {
	       ypy = double(npiy) - nq - 0.5;
	       putpixel(np,ypy,ipen);
	       putpixel(np,nq,ipen);
	     }
	     else
	       putpixel(np,nq,33);
	   }
	 } while (r <= maxiter && k<=maxiter);
       }
       if(kbhit()) break;
     }
//
//   Clean-up and end.
//
     getch();
     closegraph();
}