n-Body Solar System Software by James Pate Williams, Jr.

I have solved two planetary models of our solar system. The first was by approximating the Taylor series expansion for a n-body gravitational system. The approximate equations are given below:

n-body-problemDownload

The second computation used a fifth order Runge-Kutta method to solve the system of 3n second order non-linear ordinary differential equations. See the code below:

using System;
using System.Collections.Generic;
using System.ComponentModel;

namespace nBodyProblem
{
    public struct Vector3
    {
        public double x, y, z;
    }

    class nBodyTaylorSeries
    {
        private double G = 6.67384e-11;    // gravitational constant
        private int n;                     // number of masses
        private List<double> Mass;         // masses
        private List<Vector3> X0;          // initial positions
        private List<Vector3> V0;          // initial velocities

        public nBodyTaylorSeries(
            int n,
            List<double> Mass,
            List<Vector3> X0,
            List<Vector3> V0)
        {
            this.n = n;
            this.Mass = Mass;
            this.X0 = X0;
            this.V0 = V0;
        }

        private double Distance(Vector3 u, Vector3 v)
        {
            // Euclidean distance
            double uvxs = Math.Pow(u.x - v.x, 2);
            double uvys = Math.Pow(u.y - v.y, 2);
            double uvzs = Math.Pow(u.z - v.z, 2);

            return Math.Sqrt(uvxs + uvys + uvzs);
        }

        private Vector3 D2(
            int i,
            List<Vector3> X)
        {
            // second derivative at t0

            double sumX = 0, sumY = 0, sumZ = 0;
            Vector3 d2 = new Vector3();

            for (int k = 0; k < n; k++)
            {
                if (k != i)
                {
                    double d = Math.Pow(Distance(X[k], X[i]), 3);
                    double m = Mass[k];

                    sumX += m * (X[k].x - X[i].x) / d;
                    sumY += m * (X[k].y - X[i].y) / d;
                    sumZ += m * (X[k].z - X[i].z) / d;
                }
            }

            d2.x = G * sumX;
            d2.y = G * sumY;
            d2.z = G * sumZ;
            return d2;
        }

        private Vector3 D3(
            int i,
            List<Vector3> X,
            List<Vector3> V)
        {
            // third derivative at t0

            double sumX1 = 0, sumY1 = 0, sumZ1 = 0;
            double sumX2 = 0, sumY2 = 0, sumZ2 = 0;
            Vector3 d3 = new Vector3();

            for (int k = 0; k < n; k++)
            {
                if (k != i)
                {
                    double d = Math.Pow(Distance(X[k], X[i]), 3);
                    double m = Mass[k];

                    sumX1 += m * (V[k].x - V[i].x) / d;
                    sumY1 += m * (V[k].y - V[i].y) / d;
                    sumZ1 += m * (V[k].z - V[i].z) / d;
                }
            }

            for (int k = 0; k < n; k++)
            {
                if (k != i)
                {
                    double d = Math.Pow(Distance(X[k], X[i]), 5);
                    double m = Mass[k];

                    sumX2 += m * (X[k].x - X[i].x) * (V[k].x - V[i].x) / d;
                    sumY2 += m * (X[k].y - X[i].y) * (V[k].y - V[i].y) / d;
                    sumZ2 += m * (X[k].z - X[i].z) * (V[k].z - V[i].z) / d;
                }
            }

            d3.x = G * (sumX1 - 1.5 * sumX2);
            d3.y = G * (sumY1 - 1.5 * sumY2);
            d3.z = G * (sumZ1 - 1.5 * sumZ2);
            return d3;
        }

        private Vector3 D4(
            int i,
            List<Vector3> X,
            List<Vector3> V,
            List<Vector3> A)
        {
            // fourth derivative at t0

            double sumX1 = 0, sumY1 = 0, sumZ1 = 0;
            double sumX2 = 0, sumY2 = 0, sumZ2 = 0;
            double sumX3 = 0, sumY3 = 0, sumZ3 = 0;
            double sumX4 = 0, sumY4 = 0, sumZ4 = 0;
            Vector3 d4 = new Vector3();

            for (int k = 0; k < n; k++)
            {
                if (k != i)
                {
                    double d = Math.Pow(Distance(X[k], X[i]), 3);
                    double m = Mass[k];

                    sumX1 += m * (A[k].x - A[i].x) / d;
                    sumY1 += m * (A[k].y - A[i].y) / d;
                    sumZ1 += m * (A[k].z - A[i].z) / d;
                }
            }

            for (int k = 0; k < n; k++)
            {
                if (k != i)
                {
                    double d = Math.Pow(Distance(X[k], X[i]), 5);
                    double m = Mass[k];

                    sumX2 += m * (X[k].x - X[i].x) * (A[k].x - A[i].x) / d;
                    sumY2 += m * (X[k].y - X[i].y) * (A[k].y - A[i].y) / d;
                    sumZ2 += m * (X[k].z - X[i].z) * (A[k].z - A[i].z) / d;
                }
            }

            for (int k = 0; k < n; k++)
            {
                if (k != i)
                {
                    double d = Math.Pow(Distance(X[k], X[i]), 5);
                    double m = Mass[k];

                    sumX3 += m * Math.Pow(V[k].x - V[i].x, 2) / d;
                    sumY3 += m * Math.Pow(V[k].y - V[i].y, 2) / d;
                    sumZ3 += m * Math.Pow(V[k].z - V[i].z, 2) / d;
                }
            }

            for (int k = 0; k < n; k++)
            {
                if (k != i)
                {
                    double d = Math.Pow(Distance(X[k], X[i]), 7);
                    double m = Mass[k];

                    sumX4 += m * (X[k].x - X[i].x) * Math.Pow(V[k].x - V[i].x, 2) / d;
                    sumY4 += m * (X[k].y - X[i].y) * Math.Pow(V[k].y - V[i].y, 2) / d;
                    sumZ4 += m * (X[k].z - X[i].z) * Math.Pow(V[k].z - V[i].z, 2) / d;
                }
            }

            d4.x = G * (sumX1 - 1.5 * sumX2 - 3 * sumX3 + 15 * sumX4 / 4);
            d4.y = G * (sumY1 - 1.5 * sumY2 - 3 * sumY3 + 15 * sumY4 / 4);
            d4.z = G * (sumZ1 - 1.5 * sumZ2 - 3 * sumZ3 + 15 * sumZ4 / 4);
            return d4;
        }

        public Vector3 Position(
            int day0,
            int day1,
            int i,
            BackgroundWorker worker,
            DoWorkEventArgs e)
        {
            double t0 = day0 * 24 * 3600;
            double t1 = day1 * 24 * 3600;
            double delta = 60, t = t0;
            List<Vector3> Q0 = new List<Vector3>();
            List<Vector3> R0 = new List<Vector3>();

            for (int j = 0; j < n; j++)
            {
                Q0.Add(X0[j]);
                R0.Add(V0[j]);
            }

            while (t <= t1)
            {
                if (worker.CancellationPending)
                {
                    e.Cancel = true;
                    t = t1 + delta;
                }

                else
                {
                    List<Vector3> Q1 = new List<Vector3>();
                    List<Vector3> R1 = new List<Vector3>();
                    List<Vector3> d2 = new List<Vector3>();
                    List<Vector3> d3 = new List<Vector3>();
                    List<Vector3> d4 = new List<Vector3>();

                    for (int j = 0; j < n; j++)
                    {
                        // update derivatives
                        d2.Add(D2(j, Q0));
                        d3.Add(D3(j, Q0, R0));
                        Q1.Add(new Vector3());
                        R1.Add(new Vector3());
                    }

                    for (int j = 0; j < n; j++)
                        d4.Add(D4(i, Q0, R0, d2));

                    for (int j = 0; j < n; j++)
                    {
                        double m = Mass[j];
                        Vector3 Qj = new Vector3();
                        // update positions
                        Qj.x = Q0[j].x + delta * (R0[j].x + (delta * (d2[j].x / 2 + delta * (d3[j].x / 6 + delta * d4[j].x / 24)) / m));
                        Qj.y = Q0[j].y + delta * (R0[j].y + (delta * (d2[j].y / 2 + delta * (d3[j].y / 6 + delta * d4[j].y / 24)) / m));
                        Qj.z = Q0[j].z + delta * (R0[j].z + (delta * (d2[j].z / 2 + delta * (d3[j].z / 6 + delta * d4[j].z / 24)) / m));
                        Q1[j] = Qj;
                    }

                    for (int j = 0; j < n; j++)
                    {
                        double m = Mass[j];
                        Vector3 Rj = new Vector3();
                        // update velocities
                        Rj.x = R0[j].x + delta * (d2[j].x + delta * (d3[j].x / 2 + delta * d4[j].x / 6)) / m;
                        Rj.y = R0[j].y + delta * (d2[j].y + delta * (d3[j].y / 2 + delta * d4[j].y / 6)) / m;
                        Rj.z = R0[j].z + delta * (d2[j].z + delta * (d3[j].z / 2 + delta * d4[j].z / 6)) / m;
                        R1[j] = Rj;
                    }

                    for (int j = 0; j < n; j++)
                    {
                        Q0[j] = Q1[j];
                        R0[j] = R1[j];
                    }

                    t += delta;

                    int percentProgress = (int)(100 * t / t1);

                    worker.ReportProgress(percentProgress);
                }
            }

            // update positions and velocities
            X0 = new List<Vector3>();
            V0 = new List<Vector3>();

            for (int j = 0; j < n; j++)
            {
                X0.Add(Q0[j]);
                V0.Add(R0[j]);
            }

            return Q0[i];
        }
    }
}


using System;
using System.Collections.Generic;

namespace nBodyProblem
{
    public struct Vector3
    {
        public double x, y, z;
    }

    class RK3N
    {
        // function rk3n from "A Numerical Library in C for
        // Scientists and Engineers" by J.T. Lau PhD p. 465
        private void rk3n(ref double x, double a, double b,
            double[] y, double[] ya, double[] z, double[] za,
            Func<int, int, double, double[], double> fxyj,
            double[] e, double[] d, bool fi, int n)
        {
            bool first, last=false, reject, test, ta, tb;
            int j, jj;
            double xl, h, hmin, ind, hl=0, absh, fhm, discry, discrz, toly, tolz, mu, mu1=0, fhy, fhz;

            double[] yl = new double[n + 1];
            double[] zl = new double[n + 1];
            double[] k0 = new double[n + 1];
            double[] k1 = new double[n + 1];
            double[] k2 = new double[n + 1];
            double[] k3 = new double[n + 1];
            double[] k4 = new double[n + 1];
            double[] k5 = new double[n + 1];
            double[] ee = new double[4 * n + 1];

            if (fi)
            {
                d[3] = a;
                for (jj = 1; jj <= n; jj++)
                {
                    d[jj + 3] = ya[jj];
                    d[n + jj + 3] = za[jj];
                }
            }
            d[1] = 0.0;
            xl = d[3];
            for (jj = 1; jj <= n; jj++)
            {
                yl[jj] = d[jj + 3];
                zl[jj] = d[n + jj + 3];
            }
            if (fi) d[2] = b - d[3];
            absh = h = Math.Abs(d[2]);
            if (b - xl < 0.0) h = -h;
            ind = Math.Abs(b - xl);
            hmin = ind * e[1] + e[2];
            for (jj = 2; jj <= 2 * n; jj++)
            {
                hl = ind * e[2 * jj - 1] + e[2 * jj];
                if (hl < hmin) hmin = hl;
            }
            for (jj = 1; jj <= 4 * n; jj++) ee[jj] = e[jj] / ind;
            first = reject = true;
            test = true;
            if (fi)
            {
                last = true;
                test = false;
            }
            while (true)
            {
                if (test)
                {
                    absh = Math.Abs(h);
                    if (absh < hmin)
                    {
                        h = (h > 0.0) ? hmin : -hmin;
                        absh = hmin;
                    }
                    ta = (h >= b - xl);
                    tb = (h >= 0.0);
                    if ((ta && tb) || (!(ta || tb)))
                    {
                        d[2] = h;
                        last = true;
                        h = b - xl;
                        absh = Math.Abs(h);
                    }
                    else
                        last = false;
                }
                test = true;
                if (reject)
                {
                    x = xl;
                    for (jj = 1; jj <= n; jj++) y[jj] = yl[jj];
                    for (j = 1; j <= n; j++) k0[j] = fxyj(n, j, x, y) * h;
                }
                else
                {
                    fhy = h / hl;
                    for (jj = 1; jj <= n; jj++) k0[jj] = k5[jj] * fhy;
                }
                x = xl + 0.276393202250021 * h;
                for (jj = 1; jj <= n; jj++)
                    y[jj] = yl[jj] + (zl[jj] * 0.276393202250021 +
                                    k0[jj] * 0.038196601125011) * h;
                for (j = 1; j <= n; j++) k1[j] = fxyj(n, j, x, y) * h;
                x = xl + 0.723606797749979 * h;
                for (jj = 1; jj <= n; jj++)
                    y[jj] = yl[jj] + (zl[jj] * 0.723606797749979 +
                                    k1[jj] * 0.261803398874989) * h;
                for (j = 1; j <= n; j++) k2[j] = fxyj(n, j, x, y) * h;
                x = xl + h * 0.5;
                for (jj = 1; jj <= n; jj++)
                    y[jj] = yl[jj] + (zl[jj] * 0.5 + k0[jj] * 0.046875 + k1[jj] *
                            0.079824155839840 - k2[jj] * 0.001699155839840) * h;
                for (j = 1; j <= n; j++) k4[j] = fxyj(n, j, x, y) * h;
                x = (last ? b : xl + h);
                for (jj = 1; jj <= n; jj++)
                    y[jj] = yl[jj] + (zl[jj] + k0[jj] * 0.309016994374947 +
                                    k2[jj] * 0.190983005625053) * h;
                for (j = 1; j <= n; j++) k3[j] = fxyj(n, j, x, y) * h;
                for (jj = 1; jj <= n; jj++)
                    y[jj] = yl[jj] + (zl[jj] + k0[jj] * 0.083333333333333 + k1[jj] *
                            0.301502832395825 + k2[jj] * 0.115163834270842) * h;
                for (j = 1; j <= n; j++) k5[j] = fxyj(n, j, x, y) * h;
                reject =false;
                fhm = 0.0;
                for (jj = 1; jj <= n; jj++)
                {
                    discry = Math.Abs((-k0[jj] * 0.5 + k1[jj] * 1.809016994374947 +
                                k2[jj] * 0.690983005625053 - k4[jj] * 2.0) * h);
                    discrz = Math.Abs((k0[jj] - k3[jj]) * 2.0 - (k1[jj] + k2[jj]) * 10.0 +
                                k4[jj] * 16.0 + k5[jj] * 4.0);
                    toly = absh * (Math.Abs(zl[jj]) * ee[2 * jj - 1] + ee[2 * jj]);
                    tolz = Math.Abs(k0[jj]) * ee[2 * (jj + n) - 1] + absh * ee[2 * (jj + n)];
                    reject = ((discry > toly) || (discrz > tolz) || reject);
                    fhy = discry / toly;
                    fhz = discrz / tolz;
                    if (fhz > fhy) fhy = fhz;
                    if (fhy > fhm) fhm = fhy;
                }
                mu = 1.0 / (1.0 + fhm) + 0.45;
                if (reject)
                {
                    if (absh <= hmin)
                    {
                        d[1] += 1.0;
                        for (jj = 1; jj <= n; jj++)
                        {
                            y[jj] = yl[jj];
                            z[jj] = zl[jj];
                        }
                        first = true;
                        if (b == x) break;
                        xl = x;
                        for (jj = 1; jj <= n; jj++)
                        {
                            yl[jj] = y[jj];
                            zl[jj] = z[jj];
                        }
                    }
                    else
                        h *= mu;
                }
                else
                {
                    if (first)
                    {
                        first = false;
                        hl = h;
                        h *= mu;
                    }
                    else
                    {
                        fhy = mu * h / hl + mu - mu1;
                        hl = h;
                        h *= fhy;
                    }
                    mu1 = mu;
                    for (jj = 1; jj <= n; jj++)
                        z[jj] = zl[jj] + (k0[jj] + k3[jj]) * 0.083333333333333 +
                            (k1[jj] + k2[jj]) * 0.416666666666667;
                    if (b == x) break;
                    xl = x;
                    for (jj = 1; jj <= n; jj++)
                    {
                        yl[jj] = y[jj];
                        zl[jj] = z[jj];
                    }
                }
            }
            if (!last) d[2] = h;
            d[3] = x;
            for (jj = 1; jj <= n; jj++)
            {
                d[jj + 3] = y[jj];
                d[n + jj + 3] = z[jj];
            }
        }

        // Euclidean distance
        private double distance(
            double x1,
            double y1,
            double z1,
            double x2,
            double y2,
            double z2)
        {
            double xs = Math.Pow(x1 - x2, 2);
            double ys = Math.Pow(y1 - y2, 2);
            double zs = Math.Pow(z1 - z2, 2);

            return Math.Sqrt(xs + ys + zs);
        }

        // n-body gravitational force
        private double ftxj(int n, int j, double t, double[] x)
        {
            double G = 6.67384e-11; // gravitational constant
            double sum = 0;
            int j3 = j / 3 + 1;
            int m3 = j % 3;

            // remember we are now working with 3n equations
            for (int k = 1; k <= n / 3; k++)
            {
                int k3 = 3 * (k - 1) + 1;

                double d = distance(
                    x[k3],
                    x[k3 + 1],
                    x[k3 + 2],
                    x[j3],
                    x[j3 + 1],
                    x[j3 + 2]);
                double denom = d * d * d;

                if (denom != 0)
                    sum += (x[k3 + m3] - x[j3 + m3]) / denom;
            }

            return G * sum;
        }

        public Vector3 nBody(
            int day0,
            int day1,
            int i,
            int n,
            List<double> Mass,
            ref List<Vector3> X0,
            ref List<Vector3> V0)
        {
            int n3 = 3 * n;                 // 3n equations n (x, y, z)
            double a = day0 * 24 * 3600;    // convert to seconds
            double b = day1 * 24 * 3600;    // convert to seconds
            double x = a;
            double[] y = new double[n3 + 1];
            double[] ya = new double[n3 + 1];
            double[] z = new double[n3 + 1];
            double[] za = new double[n3 + 1];
            double[] e = new double[4 * n3 + 1];
            double[] d = new double[2 * n3 + 4];

            // initialize error constraints
            for (int j = 0; j <= 4 * n3; j++)
                e[j] = 1.0e-12;

            // create a system of 3n equations
            for (int j = 0; j < n; j++)
            {
                int j3 = 3 * j;

                ya[j3 + 1] = X0[j].x;
                ya[j3 + 2] = X0[j].y;
                ya[j3 + 3] = X0[j].z;
                za[j3 + 1] = V0[j].x;
                za[j3 + 2] = V0[j].y;
                za[j3 + 3] = V0[j].z;
            }

            // solve the system of 3n ordinary differential equations
            rk3n(ref x, a, b, y, ya, z, za, ftxj, e, d, true, n3);
            // we return the position of the desired ith body
            Vector3 result = new Vector3();

            int i3 = 3 * i;

            result.x = y[i3 + 1];
            result.y = y[i3 + 2];
            result.z = y[i3 + 3];

            // update the positions and velocities
            // wipe the slates clean
            X0 = new List<Vector3>();
            V0 = new List<Vector3>();

            for (int j = 0; j < n; j++)
            {
                int j3 = 3 * j;
                Vector3 X = new Vector3();  // new position
                Vector3 V = new Vector3();  // new velocity

                X.x = y[j3 + 1];    // 3n position components
                X.y = y[j3 + 2];
                X.z = y[j3 + 3];
                V.x = z[j3 + 1];    // 3n velocity components
                V.y = z[j3 + 2];
                V.z = z[j3 + 3];
                // update vectors
                X0.Add(X);
                V0.Add(V);
            }

            return result;
        }
    }
}

Results will be added at a later date.

Unknown's avatar

Author: jamespatewilliamsjr

My whole legal name is James Pate Williams, Jr. I was born in LaGrange, Georgia approximately 70 years ago. I barely graduated from LaGrange High School with low marks in June 1971. Later in June 1979, I graduated from LaGrange College with a Bachelor of Arts in Chemistry with a little over a 3 out 4 Grade Point Average (GPA). In the Spring Quarter of 1978, I taught myself how to program a Texas Instruments desktop programmable calculator and in the Summer Quarter of 1978 I taught myself Dayton BASIC (Beginner's All-purpose Symbolic Instruction Code) on LaGrange College's Data General Eclipse minicomputer. I took courses in BASIC in the Fall Quarter of 1978 and FORTRAN IV (Formula Translator IV) in the Winter Quarter of 1979. Professor Kenneth Cooper, a genius poly-scientist taught me a course in the Intel 8085 microprocessor architecture and assembly and machine language. We would hand assemble our programs and insert the resulting machine code into our crude wooden box computer which was designed and built by Professor Cooper. From 1990 to 1994 I earned a Bachelor of Science in Computer Science from LaGrange College. I had a 4 out of 4 GPA in the period 1990 to 1994. I took courses in C, COBOL, and Pascal during my BS work. After graduating from LaGrange College a second time in May 1994, I taught myself C++. In December 1995, I started using the Internet and taught myself client-server programming. I created a website in 1997 which had C and C# implementations of algorithms from the "Handbook of Applied Cryptography" by Alfred J. Menezes, et. al., and some other cryptography and number theory textbooks and treatises.

Leave a comment