/* ***********************************************************************
   *									 *
   *	                Quaternion Fractal Generator			 *
   *									 *
   *		     	    f(q) = q * q + lambda			 *
   *									 *
   *		      Iteration of the Forward Equation			 *
   *			     in Quaternion Space			 *
   *									 *
   *				Program by				 *
   *			   Christopher D. Watkins			 *
   *									 *
   *			     'C' conversion by				 *
   *                            Larry Sharp				 *
   *									 *
   ***********************************************************************

*/

#include "stdio.h"
#include "stdlib.h"
#include "dos.h"
#include "conio.h"
#include "math.h"
#include "string.h"
#include "math.inc"
#include "graph.inc"
#include "render.inc"
#include "quater.inc"

FDA q, l;
int DivergenceBoundary;
int MaximumIterations;
float XLeft, XRight;
float YBottom, YTop;
float ZBack, ZFront;

/* ***********************************************************************
   *									 *
   *		      Constants for Quaternion Generator		 *
   *									 *
   ***********************************************************************

   QSetCons - set up constants for quaternion generator
*/


char ObjF[]="DRAGON";

void QSetCons()
{
  q[0]=0.0;
  q[1]=0.0;
  q[2]=0.0;
  q[3]=0.0;
  l[0]=0.2809;
  l[1]=0.53;
  l[2]=0.0;
  l[3]=0.0;
  DivergenceBoundary=4;
  MaximumIterations=100;
  XLeft=-1.3;
  XRight=1.3;
  YBottom=-1.3;
  YTop=1.3;
  ZBack=0.0;
  ZFront=1.3;
}

/*
char ObjF[]="REVSOLID";

void QSetCons()
{
  q[0]=0.0;
  q[1]=0.0;
  q[2]=0.0;
  q[3]=0.0;
  l[0]=-1.0;
  l[1]=0.0;
  l[2]=0.0;
  l[3]=0.0;
  DivergenceBoundary=4;
  MaximumIterations=100;
  XLeft=-1.8;
  XRight=1.8;
  YBottom=-1.8;
  YTop=1.8;
  ZBack=0.0;
  ZFront=1.8;
}
*/

float XScale, YScale, ZScale;
float XOffset, YOffset;

void Parameters()
{
  int tmp;

  textcolor(YELLOW);
  textbackground(BLUE);
  clrscr();
  printf("Quaternion Generator\n\n");
  printf("Equation : f(q)=lambda+q*q\n\n");
  printf("lambda   : quaternion constant for a particular fractal curve\n\n");
  printf("q        : quaternion variable seeded as 0+0i+0j+0k\n\n");
  printf("The fractal surface is realized by examination of the\n");
  printf("  divergence properties of the equation in quaternion space\n\n");
  ClearHeightBuffer();
  QSetCons();
  XScale=(float) Res/(XRight-XLeft);
  XOffset=-XScale*XLeft;
  YScale=(float) Res/(YBottom-YTop);
  YOffset=-YScale*YTop;
  tmp=Res*Scaling;
  ZScale=(float) tmp/(ZFront-ZBack);
}

void Func(FDA Q, FDA NQ)
{
  QSquare(Q, Q);
  QAdd(Q, l, NQ);
}

void Iterate(FDA p, int *NumberOfIterations)
{
  FDA nq;
  float tmp;

  Func(p, q);
  *NumberOfIterations=0;
  do
  {
    Func(q, nq);
    ++*NumberOfIterations;
    q[0]=nq[0];
    q[1]=nq[1];
    q[2]=nq[2];
    q[3]=nq[3];
    tmp=SqrFP(nq[0])+SqrFP(nq[1]);
  }
  while((SqrFP(nq[0])+SqrFP(nq[1])<DivergenceBoundary) && (*NumberOfIterations!=MaximumIterations));
}

Byte MaxDepth;
FDA p;
float X, Y, Z;
int NumOfIters;

void UpdateHeight(float x, float y, float z)
{
  int PlotX, PlotY, PlotZ;

  PlotX=Round(XScale*x+XOffset);
  PlotY=Round(YScale*y+YOffset);
  PlotZ=Round(ZScale*z);
  Height[PlotX][PlotY]=PlotZ;
}

void CalculateMaxDepthOfRecursion()
{
  float k;

  k=(float) Res*0.5;
  MaxDepth=0;
  while(k>=1.0)
  {
    k*=0.5;
    ++MaxDepth;
  }
  printf("Recursion Depth based on Resolution is %d\n", MaxDepth);
}

void RecursiveCalc(float Zrec, float Delta, int Depth)
{
  Z=Zrec;
  if(Depth!=MaxDepth)
  {
    p[0]=X;
    p[1]=Y;
    p[2]=Z;
    p[3]=0.0;
    Iterate(p, &NumOfIters);
    if(NumOfIters<MaximumIterations)
      Z-=Delta;
    else
      Z+=Delta;
    RecursiveCalc(Z, Delta*0.5, Depth+1);
  }
}

void ScanSpace()
{

  float dx, dy, dz;
  int Slice;
  float Zseed, DZseed;
  int Done;

  CalculateMaxDepthOfRecursion();
  dx=(XRight-XLeft)/(float) Res;
  dy=(YBottom-YTop)/(float) Res;
  dz=(ZFront-ZBack);
  Zseed=dz*0.5;
  DZseed=Zseed*0.5;
  Slice=0;
  X=XLeft;
  do
  {
    Y=YTop;
    do
    {
      Z=ZBack;
      p[0]=X;
      p[1]=Y;
      p[2]=Z;
      p[3]=0.0;
      Iterate(p, &NumOfIters);
      if(NumOfIters<MaximumIterations)
	goto Done;
      Z=ZFront;
      p[0]=X;
      p[1]=Y;
      p[2]=Z;
      p[3]=0.0;
      Iterate(p, &NumOfIters);
      if(!(NumOfIters<MaximumIterations))
	goto Done;
      RecursiveCalc(Zseed, DZseed, 0);
Done:
      UpdateHeight(X, Y, Z);
      Y+=dy;
    }
    while(!(Y<YBottom));
    gotoxy(1, 20);
    printf("Slice # %d\n", Slice);
    ++Slice;
    X+=dx;
  }
  while(!(X>XRight));
}

/* ***********************************************************************
   *									 *
   *				Main Program				 *
   *									 *
   ***********************************************************************
*/

void main()
{
  Res=MaxRes;
  MaxHeight=Res;
  HeightBufferScalingFactor();
  Parameters();
  ScanSpace();
  strcpy(ObjectFile, ObjF);
  SaveHeightBuffer(ObjectFile);
}