//	Spludlow Software
//	Copyright © Samuel P. Ludlow 2020 All Rights Reserved
//	Distributed under the terms of the GNU General Public License version 3
//	Distributed WITHOUT ANY WARRANTY; without implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE
//	https://www.spludlow.co.uk/LICENCE.TXT
//	The Spludlow logo is a registered trademark of Samuel P. Ludlow and may not be used without permission
//	v1.14

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using System.Numerics;

namespace Spludlow.Drawing
{
	public class FractalMaths
	{
		//	Passing a null info means return the default extra paramter values

		//	Each maths method goes through each line and uses the Spludlow.Drawing.Fractals.OnTheLine() method to determine if this host/thread should render the line

		private static double[] ParseParamters(string text)
		{
			string[] words = Spludlow.Text.Split(text, '/', true, false);
			double[] result = new double[words.Length];
			for (int index = 0; index < words.Length; ++index)
				result[index] = Double.Parse(words[index]);
			return result;
		}

		public static string Mandelbrot(Spludlow.Drawing.Fractals.RenderInfo info)
		{
			if (info == null)
				return "4.0";

			double escapeRadius = Double.Parse(info.Fractal.Paramters);

			double x1 = info.Fractal.X - info.Fractal.Width / 2.0;
			double y1 = info.Fractal.Y - info.Fractal.Height / 2.0;
			double xd = info.Fractal.Width / (double)info.Fractal.PixelWidth;
			double yd = info.Fractal.Height / (double)info.Fractal.PixelHeight;

			for (int yPixel = 0; yPixel < info.Fractal.PixelHeight; yPixel++)
			{
				if (Spludlow.Drawing.Fractals.OnTheLine(yPixel, info) == false)
					continue;

				double y = y1 + yPixel * yd;

				for (int xPixel = 0; xPixel < info.Fractal.PixelWidth; xPixel++)
				{
					double x = x1 + xPixel * xd;

					int iteration = 0;
					double r1 = 0;
					double i1 = 0;
					double r1pow2 = 0;
					double i1pow2 = 0;

					double rpow = 0;
					double rlastpow = 0;

					while ((iteration < info.Fractal.Iterations) && (rpow < escapeRadius))
					{
						r1pow2 = r1 * r1;
						i1pow2 = i1 * i1;
						i1 = 2 * i1 * r1 + y;
						r1 = r1pow2 - i1pow2 + x;
						rlastpow = rpow;
						rpow = r1pow2 + i1pow2;
						iteration++;
					}

					if (iteration < info.Fractal.Iterations)
						Spludlow.Drawing.Fractals.Plot(xPixel, yPixel, iteration, info);
				}
			}

			return null;
		}

		public static string Mandelbrot2(Spludlow.Drawing.Fractals.RenderInfo info)
		{
			if (info == null)
				return "4.0";

			decimal escapeRadius = Decimal.Parse(info.Fractal.Paramters);

			decimal x1 = (decimal)info.Fractal.X - (decimal)info.Fractal.Width / 2.0M;
			decimal y1 = (decimal)info.Fractal.Y - (decimal)info.Fractal.Height / 2.0M;
			decimal xd = (decimal)info.Fractal.Width / info.Fractal.PixelWidth;
			decimal yd = (decimal)info.Fractal.Height / info.Fractal.PixelHeight;

			for (int yPixel = 0; yPixel < info.Fractal.PixelHeight; yPixel++)
			{
				if (Spludlow.Drawing.Fractals.OnTheLine(yPixel, info) == false)
					continue;

				decimal y = y1 + yPixel * yd;

				for (int xPixel = 0; xPixel < info.Fractal.PixelWidth; xPixel++)
				{
					decimal x = x1 + xPixel * xd;

					int iteration = 0;

					decimal z_r = 0;
					decimal z_i = 0;
					decimal zrsqr = 0;
					decimal zisqr = 0;
					while ((iteration < info.Fractal.Iterations) && (zrsqr + zisqr <= escapeRadius))
					{
						z_i = z_r * z_i;
						z_i += z_i;
						z_i += y;
						z_r = zrsqr - zisqr + x;
						zrsqr = z_r * z_r;
						zisqr = z_i * z_i;

						++iteration;
					}

					if (iteration < info.Fractal.Iterations)
						Spludlow.Drawing.Fractals.Plot(xPixel, yPixel, iteration, info);
				}
			}

			return null;
		}

		public static string Julia(Spludlow.Drawing.Fractals.RenderInfo info)
		{
			if (info == null)
				return "-0.7 / 0.3";

			double[] paramters = ParseParamters(info.Fractal.Paramters);

			if (paramters.Length != 2)
				throw new ApplicationException("FractalMaths, Julia; Expected 2 paramters:\t" + info.Fractal.Paramters);

			double xc = paramters[0];
			double yc = paramters[1];

			double x1 = info.Fractal.X - info.Fractal.Width / 2.0;
			double y1 = info.Fractal.Y - info.Fractal.Height / 2.0;
			double xd = info.Fractal.Width / (double)info.Fractal.PixelWidth;
			double yd = info.Fractal.Height / (double)info.Fractal.PixelHeight;

			for (int yPixel = 0; yPixel < info.Fractal.PixelHeight; yPixel++)
			{
				if (Spludlow.Drawing.Fractals.OnTheLine(yPixel, info) == false)
					continue;

				double y = y1 + yPixel * yd;

				for (int xPixel = 0; xPixel < info.Fractal.PixelWidth; xPixel++)
				{
					double x = x1 + xPixel * xd;

					int iteration = 0;

					double r1 = x;
					double i1 = y;
					double r1pow2 = x * x;
					double i1pow2 = y * y;

					double rpow = 0;
					double rlastpow = 0;

					while ((iteration < info.Fractal.Iterations) && (rpow < 4))		//	4
					{
						r1pow2 = r1 * r1;
						i1pow2 = i1 * i1;
						i1 = 2 * i1 * r1 + yc;
						r1 = r1pow2 - i1pow2 + xc;
						rlastpow = rpow;
						rpow = r1pow2 + i1pow2;
						iteration++;
					}

					if (iteration < info.Fractal.Iterations)
						Spludlow.Drawing.Fractals.Plot(xPixel, yPixel, iteration, info);
				}
			}

			return null;
		}

		public static string Newton(Spludlow.Drawing.Fractals.RenderInfo info)
		{
			if (info == null)
				return "1.0";

			double a = Double.Parse(info.Fractal.Paramters);
			
			double x1 = info.Fractal.X - info.Fractal.Width / 2.0;
			double y1 = info.Fractal.Y - info.Fractal.Height / 2.0;
			double xd = info.Fractal.Width / (double)info.Fractal.PixelWidth;
			double yd = info.Fractal.Height / (double)info.Fractal.PixelHeight;

			for (int yPixel = 0; yPixel < info.Fractal.PixelHeight; yPixel++)
			{
				if (Spludlow.Drawing.Fractals.OnTheLine(yPixel, info) == false)
					continue;

				double y = y1 + yPixel * yd;

				for (int xPixel = 0; xPixel < info.Fractal.PixelWidth; xPixel++)
				{
					double x = x1 + xPixel * xd;

					if (x == 0 && y == 0)
						continue;

					Complex zn = new Complex(x, y);
					Complex pz = new Complex(1, 0);

					int iteration = 0;
					while ((iteration < info.Fractal.Iterations) && Spludlow.Maths.ModulusSquared(pz) > 0.00000001)
					{
						pz = Spludlow.Maths.Pow(zn, 3) - 2;

						Complex pzd = 3 * Spludlow.Maths.Pow(zn, 2);
						zn = (zn - a * pz / pzd);
						iteration++;
					}

					if (iteration < info.Fractal.Iterations)
						Spludlow.Drawing.Fractals.Plot(xPixel, yPixel, iteration, info);
				}
			}

			return null;
		}

		public static string NewtonRoot(Spludlow.Drawing.Fractals.RenderInfo info)
		{
			if (info == null)
				return "5";

			int rootCount = Int32.Parse(info.Fractal.Paramters);

			double x1 = info.Fractal.X - info.Fractal.Width / 2.0;
			double y1 = info.Fractal.Y - info.Fractal.Height / 2.0;
			double xd = info.Fractal.Width / (double)info.Fractal.PixelWidth;
			double yd = info.Fractal.Height / (double)info.Fractal.PixelHeight;

			Complex[] roots = new Complex[rootCount];
			for (int index = 0; index < rootCount; ++index)
			{
				double r = Math.Cos(2 * index * Math.PI / rootCount);
				double i = Math.Sin(2 * index * Math.PI / rootCount);
				roots[index] = new Complex(r, i);
			}

			double cutoff = 0.00000000001;

			for (int yPixel = 0; yPixel < info.Fractal.PixelHeight; yPixel++)
			{
				if (Spludlow.Drawing.Fractals.OnTheLine(yPixel, info) == false)
					continue;

				double y = y1 + yPixel * yd;

				for (int xPixel = 0; xPixel < info.Fractal.PixelWidth; xPixel++)
				{
					double x = x1 + xPixel * xd;

					Complex complex = new Complex(x, y);
					Complex epsilon;

					int iteration = 0;
					do
					{
						if (++iteration > info.Fractal.Iterations)
							break;

						Complex F = Spludlow.Maths.Power(complex, new Complex(rootCount, 0)) - 1;

						Complex dFdx = rootCount * Spludlow.Maths.Power(complex, new Complex(rootCount - 1, 0));
									
						epsilon = -(F / dFdx);

						complex += epsilon;
					}
					while (Spludlow.Maths.ModulusSquared(epsilon) > cutoff);

					if (iteration < info.Fractal.Iterations)
					{
						for (int r = 0; r <= roots.GetUpperBound(0); r++)
						{
							if (Spludlow.Maths.IsCloseTo(complex, roots[r]) == true)
							{
								Spludlow.Drawing.Fractals.Plot(xPixel, yPixel, iteration, info);
								break;
							}
						}
					}
				}
			}

			return null;
		}


		public static string NewtonRhapson(Spludlow.Drawing.Fractals.RenderInfo info)
		{
			if (info == null)
				return "1.2";

			double denomMultiply = Double.Parse(info.Fractal.Paramters);

			double x1 = info.Fractal.X - info.Fractal.Width / 2.0;
			double y1 = info.Fractal.Y - info.Fractal.Height / 2.0;
			double xd = info.Fractal.Width / (double)info.Fractal.PixelWidth;
			double yd = info.Fractal.Height / (double)info.Fractal.PixelHeight;
			
			for (int yPixel = 0; yPixel < info.Fractal.PixelHeight; yPixel++)
			{
				if (Spludlow.Drawing.Fractals.OnTheLine(yPixel, info) == false)
					continue;

				double startY = y1 + yPixel * yd;

				for (int xPixel = 0; xPixel < info.Fractal.PixelWidth; xPixel++)
				{
					double y = startY;
					double x = x1 + xPixel * xd;

					double yOld = 42;
					double xOld = 42;

					int iteration;
					for (iteration = 0; iteration < info.Fractal.Iterations; ++iteration)
					{
						double Xsquare = x * x;
						double Ysquare = y * y;

						double denom = 3 * ((Xsquare - Ysquare) * (Xsquare - Ysquare) + 4 * Xsquare * Ysquare);
						if (denom == 0)
							denom = .0000001;

						x = 0.6666667 * x + (Xsquare - Ysquare) / (denom * 1);
						y = 0.6666667 * y - 2 * x * y / (denom * denomMultiply);

						if ((xOld == x) && (yOld == y))
							break;
						xOld = x;
						yOld = y;
					}

					if (iteration < info.Fractal.Iterations)
						Spludlow.Drawing.Fractals.Plot(xPixel, yPixel, iteration, info);
				}
			}

			return null;
		}


	}
}