using System;

namespace Muschel_Kreuzungen
{
    internal class MuschelMaler
    {
        private const double Ba = (double)0.0;
        private const double Bb = (double)0.01;
        private const double Bc = (double)0.02;
        private const double Da = (double)0.005;
        private const double Db = (double)0.0;
        private const double Dc = (double)0.3;
        private const double Ra = (double)0.0;
        private const double Rb = (double)0.0;
        private const double Rc = (double)0.0;
        private const double Sa = (double)0.05;
        
        public static Array[] Male(int xWidth, double tMax, double tDelta)
        {
            Random rand = new Random();
            var steps = (int)Math.Ceiling(tMax / tDelta);
            tDelta = tMax / steps;
            var aResultArray = new double[steps, xWidth];
            var bResultArray = new double[steps, xWidth];
            var cResultArray = new double[steps, xWidth];
            double aToT, bToT, cToT;
            double h = 1.0 / ((double)xWidth - 1.0);
            
            for (int i = 0; i < xWidth; i++)
            {
                aResultArray[0, i] = 1e-3*rand.NextDouble();
                bResultArray[0, i] = 1e-3*rand.NextDouble();
                cResultArray[0, i] = 1e-3*rand.NextDouble();
            }
            for (int i = 1; i < steps; i++)
            {
                for (int j = 0; j < xWidth; j++)
                {
                    double tempDouble = Convert.ToDouble(0.9 + 0.2 * rand.NextDouble());
                    double s = Ra * tempDouble;
                    double astar, asquare, aqb, aqc;

                    asquare = aResultArray[i - 1, j] * aResultArray[i - 1, j];
                    astar = (asquare + Ba) / (1 + Sa * asquare);

                    aqb = s * astar * bResultArray[i - 1, j];
                    aqc = s * astar * cResultArray[i - 1, j];

                    // linker Rand
                    if (j == 0)
                    {
                        aToT = aqb + aqc 
                        - Ra * aResultArray[i - 1, j] +
                        Da/(h*h) * (aResultArray[i - 1, j + 1] - aResultArray[i - 1, j]);
                        bToT = Bb
                        - aqb 
                        - Rb * bResultArray[i - 1, j] +
                        Db/(h*h) * (bResultArray[i - 1, j + 1] - bResultArray[i - 1, j]);
                        cToT = Bc -
                        aqc - Rc * cResultArray[i - 1, j] +
                        Dc/(h*h) * (cResultArray[i - 1, j + 1] - cResultArray[i - 1, j]);
                    }
                    // rechter Rand
                    else if (j == xWidth - 1)
                    {
                        aToT = aqb + aqc - Ra * aResultArray[i - 1, j] +
                        Da/(h*h) * (aResultArray[i - 1, j - 1] - aResultArray[i - 1, j]);
                        bToT = Bb -
                        aqb - Rb * bResultArray[i - 1, j] +
                        Db/(h*h) * (bResultArray[i - 1, j - 1] - bResultArray[i - 1, j]);
                        cToT = Bc -
                        aqc - Rc * cResultArray[i - 1, j] +
                        Dc/(h*h) * (cResultArray[i - 1, j - 1] - cResultArray[i - 1, j]);
                    }
                    else
                    {
                        // innere Punkte
                        aToT = aqb + aqc - Ra * aResultArray[i - 1, j] +
                        Da/(h*h) * (aResultArray[i - 1, j - 1] + aResultArray[i - 1, j + 1] - 2 * aResultArray[i - 1, j]);
                        bToT = Bb -
                        aqb - Rb * bResultArray[i - 1, j] +
                        Db/(h*h) * (bResultArray[i - 1, j - 1] + bResultArray[i - 1, j + 1] - 2 * bResultArray[i - 1, j]);
                        cToT = Bc -
                        aqc - Rc * cResultArray[i - 1, j] +
                        Dc/(h*h) * (cResultArray[i - 1, j - 1] + cResultArray[i - 1, j + 1] - 2 * cResultArray[i - 1, j]);
                    }
                    
                    aResultArray[i, j] = aResultArray[i - 1, j] + aToT * tDelta;
                    bResultArray[i, j] = bResultArray[i - 1, j] + bToT * tDelta;
                    cResultArray[i, j] = aResultArray[i - 1, j] + cToT * tDelta;
                    
                }
            }
            var resultArray = new Array[3];
            resultArray[0] = aResultArray;
            resultArray[1] = bResultArray;
            resultArray[2] = cResultArray;
            return resultArray;
        }
    }
}