using System;
using System.Collections.Generic;

namespace Science.Mathematics.VectorCalculus
{
    public class LevelSurface
    {
        private Function.ToLastType<double, double, double, double> f;
        private double c;
        private Domain3D U;
        public LevelSurface(Function.ToLastType<double, double, double, double> function, double constant, Domain3D domain)
        {
            f = function;
            c = constant;
            U = domain;
        }
        public double Constant
        {
            set { c = value; }
        }

        private List<double> spx = new List<double>();
        public List<double> SurfacePointX
        {
            get { return spx; }
        }
        private List<double> spy = new List<double>();
        public List<double> SurfacePointY
        {
            get { return spy; }
        }
        private List<double> spz = new List<double>();
        public List<double> SurfacePointZ
        {
            get { return spz; }
        }

        public void FindSurface(double meshSize)
        {
            int nx = (int)((U.UpperBoundOfX - U.LowerBoundOfX) / meshSize);
            int ny = (int)((U.UpperBoundOfY - U.LowerBoundOfY) / meshSize);
            int nz = (int)((U.UpperBoundOfZ - U.LowerBoundOfZ) / meshSize);
            double x, y, z;
            double v1, v2, v3, v4;
            for (int k = 0; k < nx - 1; k++)
            {
                for (int m = 0; m < ny - 1; m++)
                {
                    for (int n = 0; n < nz - 1; n++)
                    {
                        x = U.LowerBoundOfX + meshSize * k;
                        y = U.LowerBoundOfX + meshSize * m;
                        z = U.LowerBoundOfX + meshSize * n;
                        if (U.Check(x, y, z))
                        {
                            v1 = f(x, y, z) - c;
                            v2 = f(x + meshSize, y, z) - c;
                            v3 = f(x, y + meshSize, z) - c;
                            v4 = f(x, y, z + meshSize) - c;
                            if (v1 * v2 <= 0.0)
                            {
                                spx.Add(x + meshSize * Math.Abs(v1) / (Math.Abs(v1) + Math.Abs(v2)));
                                spy.Add(y);
                                spz.Add(z);
                            }
                            if (v1 * v3 <= 0.0)
                            {
                                spx.Add(x);
                                spy.Add(y + meshSize * Math.Abs(v1) / (Math.Abs(v1) + Math.Abs(v3)));
                                spz.Add(z);
                            }
                            if (v1 * v4 <= 0.0)
                            {
                                spx.Add(x);
                                spy.Add(y);
                                spz.Add(z + meshSize * Math.Abs(v1) / (Math.Abs(v1) + Math.Abs(v4)));
                            }
                        }                   
                    }
                }
            }
        }
    }
}