//
//+----------------------------------------------------------------------+
//+ Program SIEGEL.CPP                                                   +
//+ By Fausto A. A. Barbuto, BJ06@C53000.PETROBRAS.ANRJ.BR               +
//+ Rio de Janeiro, BRAZIL, on May 13, 1994.                             +
//+                                                                      +
//+ Plots a Julia set of the type z -> z^2 + lambda*z, where lambda =    +
//+ exp(2*Pi*GM*i) and GM = (5^(1/2) - 1)/2 (Golden Mean). Siegel disks  +
//+ (are they in fact Siegel disks? I'm not sure) are shown in the inner +
//+ part of the plot.                                                    +
//+                                                                      +
//+ Case with z -> z^2 + 1.001*z*exp(2*Pi*i/20) added Sunday, May 15.    +
//+ (see Fig. 30, pg. 40)                                                +
//+                                                                      +
//+                                                                      +
//+ REFERENCE: Peitgen, H.-O. & Richter, P.H.: "The Beauty of Fractals", +
//+            Springer-Verlag, 1986, pg. 30-31, pg. 77 (Map 25).        +
//+                                                                      +
//+ Needs SVGA256.BGI and SVGA256.H (from SVGABG50 package by Jordan     +
//+ Hargrave).                                                           +
//+                                                                      +
//+ Authorized version for spanky.triumf.ca site (Vancouver, BC, CANADA).+
//+----------------------------------------------------------------------+
//
#include <stdio.h>
#include <conio.h>
#include <graphics.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>
#include "Svga256.h"

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

int Vid;  //* Global variable *//

int huge DetectSVGA256()
{
  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);
  return Vid;
}

void main(void)
{
      complex lambda, i;
      double xmin, xmax, ymin, ymax, r, R1, I1, RL, IL;
      double x,y, x0,y0, GM, deltax,deltay, s1, s2, Pi=3.1415926535897932;
      register int npix, npiy, kcolor=1024, k, np, nq, ipen;
      int graphdriver=DETECT, graphmode, icase;

      clrscr();
      printf("\n  Program SIEGEL.CPP \n");
      printf("\n\n  By Fausto A. A. Barbuto, May 14, 1994");
      printf("\n  Rio de Janeiro, Federal Republic of Brazil");
      printf("\n  E-mail: BJ06@C53000.PETROBRAS.ANRJ.BR\n\n");
      printf("\n  Reference: Peitgen, H.-O., and Richter, P.H. :");
      printf("\n  'The Beauty of Fractals', Springer-Verlag, 1986");
      printf("\n  Section 2, pp. 27-32\n\n");
      printf("\n  Enter the case:\n");
      printf("\n  (1) z -> z*z +   lambda*z  [lambda=exp(2*Pi*i*GM)]\n");
      printf("                   [GM = Golden Mean = (5^(1/2) - 1.0)/2]\n");
      printf("\n  (2) z -> z*z + s*lambda*z (lambda=exp(2*Pi*i/20), s=1.001");
      printf("\n\n > ");
      scanf("%d",&icase);
      if ((icase < 1) || (icase > 2)) icase = 1;
      clrscr();

      installuserdriver("Svga256",DetectSVGA256);
      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<0) || (Vid)>4) Vid = 2;

      GM = 0.5*(sqrt(5.0) - 1.0);  //* The Golden Mean *//
      i = complex(0.0,1.0);

      if (icase == 1) {
	lambda = exp(2.0*Pi*GM*i);
	xmin = -1.05; xmax = 1.75;
	ymin = -0.75; ymax = 1.45;
      }
      else {
	lambda = 1.001*exp(0.1*Pi*i);
	xmin = -1.75; xmax = 0.75;
	ymin = -1.40; ymax = 1.15;
      }
      RL = real(lambda);
      IL = imag(lambda);
      deltax = (xmax-xmin)/(npix-1);
      deltay = (ymax-ymin)/(npiy-1);

      cleardevice();
      for (np=0; np<=npix-1; np++) {
	x0 = xmin + (double)np*deltax;
	for (nq=0; nq<=npiy-1; nq++) {
	y0 = ymin + (double)nq*deltay;
	x = x0;
	y = y0;
	k  = 0;

	 do {
	   s1 = (x+RL);
	   s2 = (y+IL);
	   R1 = x*s1 - y*s2;
	   I1 = x*s2 + y*s1;
	   r = sqrt(R1*R1 + I1*I1);
	   k++;
//*
//*        If r >= kcolor the point escapes towards infinity.
//*
	   if (r >= kcolor) {
	     ipen = 29 + k;
	     putpixel(np,nq,ipen);
	   }
//*
//*        Converging points.
//*
	   if (k == kcolor) {
//*
//*          The colour shades of the "Siegel disks" (?) are defined here.
//*          (as function of the  distance between the present point [x,y]
//*          and the "invariant" point [x=0,y=0])
//*
	     ipen = 17 + (int)(25.0*sqrt(x*x+y*y));
	     putpixel(np,nq,ipen);
	   }
//*
//*        Returns if no convergence's achieved on either escape or atraction.
//*
	   x = R1;
	   y = I1;
	 } while (r<=kcolor && k<=kcolor);
       }
       if(kbhit()) break;
     }
//*
//*  Clean-up.
//*
     getch();
     closegraph();
}