wellPath

mincurve.cs

using System;
using System.Collections.Generic;
using System.Windows.Media.Media3D;

namespace YourNamespace
{
    public static class MinimumCurvature
    {
        public static List<Point3D> MinimumCurvatureInner(double[] md, double[] inc, double[] azi)
        {
            double[][] dv = Geometry.DirectionVectorRadians(inc, azi);
            double[][] upper = new double[dv.Length - 1][];
            double[][] lower = new double[dv.Length - 1][];
            for (int i = 0; i < dv.Length - 1; i++)
            {
                upper[i] = dv[i];
                lower[i] = dv[i + 1];
            }

            double[] dogleg = Geometry.AngleBetween(upper, lower);

            double[] rf = new double[dogleg.Length];
            for (int i = 0; i < dogleg.Length; i++)
            {
                if (dogleg[i] == 0)
                    rf[i] = 1;
                else
                    rf[i] = (2 * Math.Tan(dogleg[i] / 2)) / dogleg[i];
            }

            double[] mdDiff = new double[md.Length - 1];
            for (int i = 0; i < md.Length - 1; i++)
            {
                mdDiff[i] = md[i + 1] - md[i];
            }

            double[] halfMd = new double[mdDiff.Length];
            for (int i = 0; i < mdDiff.Length; i++)
            {
                halfMd[i] = mdDiff[i] / 2;
            }

            List<Point3D> points = new List<Point3D>();

            double northing = 0;
            double easting = 0;
            double tvd = 0;

            for (int i = 0; i < halfMd.Length; i++)
            {
                northing += halfMd[i] * (upper[i][0] + lower[i][0]) * rf[i];
                easting += halfMd[i] * (upper[i][1] + lower[i][1]) * rf[i];
                tvd += halfMd[i] * (upper[i][2] + lower[i][2]) * rf[i];

                points.Add(new Point3D(easting, northing, tvd));
            }

            return points;
        }

        public static List<Point3D> MinimumCurvature(double[] md, double[] inc, double[] azi, double courseLength = 30)
        {
            try
            {
                courseLength += 0;
            }
            catch (ArgumentException)
            {
                throw new ArgumentException("courseLength must be a float");
            }

            double[] mdCopy = new double[md.Length];
            double[] incCopy = new double[inc.Length];
            double[] aziCopy = new double[azi.Length];
            Array.Copy(md, mdCopy, md.Length);
            Array.Copy(inc, incCopy, inc.Length);
            Array.Copy(azi, aziCopy, azi.Length);

            for (int i = 0; i < incCopy.Length; i++)
            {
                incCopy[i] = MathUtils.Deg2Rad(incCopy[i]);
            }

            for (int i = 0; i < aziCopy.Length; i++)
            {
                aziCopy[i] = MathUtils.Deg2Rad(aziCopy[i]);
            }

            double[] mdDiff = new double[mdCopy.Length - 1];
            for (int i = 0; i < mdDiff.Length; i++)
            {
                mdDiff[i] = mdCopy[i + 1] - mdCopy[i];
            }

            List<Point3D> points = MinimumCurvatureInner(mdCopy, incCopy, aziCopy);

            points.Insert(0, new Point3D(0, 0, 0));

            double[] dl = new double[points.Count - 1];
            for (int i = 0; i < points.Count - 1; i++)
            {
                dl[i] = MathUtils.Rad2Deg(Geometry.AngleBetween(points[i], points[i + 1]));
            }

            double[] dls = new double[dl.Length];
            for (int i = 0; i < dl.Length; i++)
            {
                dls[i] = dl[i] * (courseLength / mdDiff[i]);
            }

            points.Insert(0, new Point3D(0, 0, 0));

            return points;
        }
    }

    public static class Geometry
    {
        public static double[][] DirectionVectorRadians(double[] inc, double[] azi)
        {
            int length = Math.Min(inc.Length, azi.Length);
            double[][] dv = new double[length][];
            for (int i = 0; i < length; i++)
            {
                double radInc = MathUtils.Deg2Rad(inc[i]);
                double radAzi = MathUtils.Deg2Rad(azi[i]);
                dv[i] = new double[] { Math.Sin(radInc) * Math.Cos(radAzi), Math.Sin(radInc) * Math.Sin(radAzi), Math.Cos(radInc) };
            }
            return dv;
        }

        public static double AngleBetween(Point3D upper, Point3D lower)
        {
            double dotProduct = upper.X * lower.X + upper.Y * lower.Y + upper.Z * lower.Z;
            return Math.Acos(dotProduct);
        }
    }

    public static class MathUtils
    {
        public static double Deg2Rad(double degrees)
        {
            return degrees * (Math.PI / 180);
        }

        public static double Rad2Deg(double radians)
        {
            return radians * (180 / Math.PI);
        }
    }
}