// Fraktale Newtona
// www.algorytm.org
// Tomasz Lubinski (c) 2009
//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop
#include "Newton.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}

#define MIN(a,b) ((a<b)?a:b)
#define MAX(a,b) ((a>b)?a:b)

//---------------------------------------------------------------------------
//Describes places to render
double ratioX, ratioY;
double minX, minY, maxX, maxY;
int downX, downY;
int cSet;
//colors
TColor colors[101];
//---------------------------------------------------------------------------

//for HSV colors
void HSV2RGB(float hue, float sat, float val, float &red, float &grn, float &blu)
{
int i;
float f, p, q, t;
if(val==0) {red=0; grn=0; blu=0;}
else{
 hue/=60;
 i = (int)(hue);
 f = hue-i;
 p = val*(1-sat);
 q = val*(1-(sat*f));
 t = val*(1-(sat*(1-f)));
 if (i==0) {red=val; grn=t; blu=p;}
 else if (i==1) {red=q; grn=val; blu=p;}
 else if (i==2) {red=p; grn=val; blu=t;}
 else if (i==3) {red=p; grn=q; blu=val;}
 else if (i==4) {red=t; grn=p; blu=val;}
 else if (i==5) {red=val; grn=p; blu=q;}
}
}

//initialize array with colors used to color different levels
void initializeColors()
{
        int level;
        float r, g, b;

        for (level=0; level<100; level++)
        {
                HSV2RGB(level*2.5, 0.85, 0.8, r, g, b);
                colors[level] = (TColor)(int)(r*255) + ((int)(g*255) << 8) + ((int)(b*255) << 16);
        }
        colors[100] = (TColor)0;
}

//type for complex numbers
typedef struct complex_type
{
        double real;
        double imaginary;
} complex_t;

//calculate squared modus of given complex c
double complexModSq(complex_t c)
{
        return c.real*c.real + c.imaginary*c.imaginary;
}

//complex addition
complex_t add(complex_t a, complex_t b)
{
        complex_t result;

        result.real =  a.real + b.real;
        result.imaginary = a.imaginary + b.imaginary;

        return result;
}

//complex substraction
complex_t sub(complex_t a, complex_t b)
{
        complex_t result;

        result.real =  a.real - b.real;
        result.imaginary = a.imaginary - b.imaginary;

        return result;
}


//complex multiplication
complex_t mul(complex_t a, complex_t b)
{
        complex_t result;

        result.real =  a.real*b.real - a.imaginary*b.imaginary;
        result.imaginary = a.real*b.imaginary + a.imaginary*b.real;

        return result;
}

//complex divide
complex_t div(complex_t a, complex_t b)
{
        complex_t result;
        double x = b.real*b.real + b.imaginary*b.imaginary;

        result.real =  (a.real*b.real + a.imaginary*b.imaginary) / x;
        result.imaginary = (a.imaginary*b.real - a.real*b.imaginary) / x;

        return result;
}

//func = a[0] + a[1]*z^1 + a[2]*z^2 + ... a[n]*z^2
complex_t func(complex_t z, complex_t *a, int n)
{
        int i;
        complex_t result;

        result.real = 0;
        result.imaginary = 0;

        for (i=n; i>=0; i--)
        {
                result = add(a[i], mul(result, z));
        }

        return result;
}

complex_t results[10];
int resultsCnt;

int findResult(complex_t a)
{
        int i;
        
        for (i=0; i<resultsCnt; i++)
        {
              if (complexModSq(sub(a, results[i])) < 0.1)
              {
                 return i;
              }
        }

        results[resultsCnt] = a;
        resultsCnt++;
        return resultsCnt;
}

//value is inside set in the returned level
int levelSet(complex_t a, complex_t p, complex_t *w, complex_t *d, int n)
{
        complex_t z, z_prev;
        int iteration;

        iteration = 0;
        z = p;

        do
        {
                z_prev = z;
                z = sub(z, mul(a, div(func(z, w, n), func(z, d, n-1))));
                iteration++;
        } while (complexModSq(sub(z_prev,z)) > 0.0001 && iteration < 100);

        if (cSet == true)
        {
                if (iteration < 100)
                {
                   return 10*findResult(z);
                }
        }

        return iteration;
}

//generate fractal
void __fastcall TForm1::Button1Click(TObject *Sender)
{
        int i, j, level;
        complex_t p, w[6], d[5], a;

        minX = StrToFloat(minx->Text);
        minY = StrToFloat(miny->Text);
        maxX = StrToFloat(maxx->Text);
        maxY = StrToFloat(maxy->Text);

        for (i=0; i<6; i++)
        {
                w[i].real = StrToFloat(aArray->Cells[1][i+1]);
                w[i].imaginary = StrToFloat(aArray->Cells[2][i+1]);
        }

        for (i=1; i<6; i++)
        {
                d[i-1].real = (double)i*StrToFloat(aArray->Cells[1][i+1]);
                d[i-1].imaginary = (double)i*StrToFloat(aArray->Cells[2][i+1]);
        }

        a.real = StrToFloat(aRe->Text);
        a.imaginary = StrToFloat(aIm->Text);

        ratioX = (maxX - minX) / Fractal->Width;
        ratioY = (maxY - minY) / Fractal->Height;

        resultsCnt = 0;
        cSet = colorSet->Checked;

        for (i=0; i<Fractal->Height; i++)
        {
                p.imaginary = i*ratioY + minY;
                for (j=0; j<Fractal->Width; j++)
                {
                        p.real = j*ratioX + minX;
                        level = levelSet(a, p, w, d, 5);
                        Fractal->Canvas->Pixels[j][i] = colors[level];
                }
        }
        Fractal->Refresh();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormCreate(TObject *Sender)
{
        minx->Text = FloatToStr(-1.5);
        maxx->Text = FloatToStr(1.5);
        miny->Text = FloatToStr(-1.5);
        maxy->Text = FloatToStr(1.5);

        aRe->Text = FloatToStr(1.0);
        aIm->Text = FloatToStr(0.0);

        aArray->Cells[1][0] = "p(Re)^n";
        aArray->Cells[2][0] = "p(Im)^n";
        aArray->Cells[0][0] = "n";
        aArray->Cells[0][1] = "0";
        aArray->Cells[0][2] = "1";
        aArray->Cells[0][3] = "2";
        aArray->Cells[0][4] = "3";
        aArray->Cells[0][5] = "4";
        aArray->Cells[0][6] = "5";

        aArray->Cells[1][1] = FloatToStr(-1.0);
        aArray->Cells[1][2] = FloatToStr(0.0);
        aArray->Cells[1][3] = FloatToStr(0.0);
        aArray->Cells[1][4] = FloatToStr(1.0);
        aArray->Cells[1][5] = FloatToStr(0.0);
        aArray->Cells[1][6] = FloatToStr(0.0);

        aArray->Cells[2][1] = FloatToStr(0.0);
        aArray->Cells[2][2] = FloatToStr(0.0);
        aArray->Cells[2][3] = FloatToStr(0.0);
        aArray->Cells[2][4] = FloatToStr(0.0);
        aArray->Cells[2][5] = FloatToStr(0.0);
        aArray->Cells[2][6] = FloatToStr(0.0);

        minX = StrToFloat(minx->Text);
        minY = StrToFloat(miny->Text);
        maxX = StrToFloat(maxx->Text);
        maxY = StrToFloat(maxy->Text);

        ratioX = (maxX - minX) / Fractal->Width;
        ratioY = (maxY - minY) / Fractal->Height;

        //initialize array with colors
        initializeColors();

        //render new fractal
        Button1Click(Sender);
}

//---------------------------------------------------------------------------
void __fastcall TForm1::FractalMouseDown(TObject *Sender,
      TMouseButton Button, TShiftState Shift, int X, int Y)
{
        downX = X;
        downY = Y;

        Selection->Width = 0;
        Selection->Height = 0;
        Selection->Visible = true;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FractalMouseUp(TObject *Sender,
      TMouseButton Button, TShiftState Shift, int X, int Y)
{
        //remove selection
        Selection->Visible = false;

        //get new range to render
        minx->Text = FloatToStr(MIN(downX, X)*ratioX + minX);
        maxx->Text = FloatToStr(MAX(downX, X)*ratioX + minX);
        miny->Text = FloatToStr(MIN(downY, Y)*ratioY + minY);
        maxy->Text = FloatToStr(MAX(downY, Y)*ratioY + minY);

        minX = StrToFloat(minx->Text);
        minY = StrToFloat(miny->Text);
        maxX = StrToFloat(maxx->Text);
        maxY = StrToFloat(maxy->Text);

        //render new fractal
        Button1Click(Sender);
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FractalMouseMove(TObject *Sender,
      TShiftState Shift, int X, int Y)
{
        //if left mouse button is held then draw selection
        if (Shift.Contains(ssLeft))
        {
                Selection->Width = abs(downX - X);
                Selection->Height = abs(downY - Y);
                Selection->Left = Fractal->Left + MIN(downX, X);
                Selection->Top = Fractal->Top + MIN(downY, Y);
        }
}

//---------------------------------------------------------------------------