Schwarzschild Orbit by James Pate Williams, Jr.

The Rosette motion of the perihelion precession of a massive object orbiting a black hole (Schwarzschild solution Einstein’s general relativity field equations is illustrated by the graphs below for varying values of eccentricity of the ellipsoidal orbit 0.00, 0.01, 0.02, 0.03, 0.04, and 0.05. The angular momentum constant is B = 1, and the mass is M = 10.

using System;
using System.Drawing;
using System.Windows.Forms;

namespace SchwarzschildOrbit
{
    public partial class GraphForm : Form
    {
        private double B, B2, B4, M, M2, M3, epsilon;

        public GraphForm(double B, double M, double epsilon)
        {
            InitializeComponent();
            this.B = B;
            this.M = M;
            this.epsilon = epsilon;
            B2 = B * B;
            B4 = B2 * B2;
            M2 = M * M;
            M3 = M * M2;
            panel1.Paint += Panel1_Paint;
        }

        private double u0(double phi)
        {
            return M * (1.0 + epsilon * Math.Cos(phi)) / B2;
        }

        private double u1(double phi)
        {
            return u0(phi) + 3.0 * M3 * (1.0 + epsilon * phi * Math.Sin(phi)
                + epsilon * epsilon * (0.5 - Math.Cos(2.0 * phi) / 6.0)) / B4;
        }

        private double X(double r, double phi)
        {
            return r * Math.Cos(phi);
        }

        private double Y(double r, double phi)
        {
            return r * Math.Sin(phi);
        }

        private void Maximums(out double maxR, out double maxPhi,
            out double XMax, out double XMin, out double YMax, out double YMin)
        {
            double phi = 0.0, r = 0.0, XC = 0.0, YC = 0.0;

            maxPhi = 0.0;
            maxR = double.MinValue;
            XMax = double.MinValue;
            YMax = double.MinValue;
            XMin = double.MaxValue;
            YMin = double.MaxValue;

            while (phi <= 8.0 * Math.PI)
            {
                r = 1.0 / u1(phi);

                if (r > maxR)
                {
                    maxR = r;
                    maxPhi = phi;
                }

                XC = X(r, phi);

                if (XC > XMax)
                    XMax = XC;

                YC = Y(r, phi);

                if (YC > YMax)
                    YMax = YC;

                if (XC < XMin)
                    XMin = XC;

                if (YC < YMin)
                    YMin = YC;

                phi += 0.001;
            }
        }

        private void Minimums(out double minR, out double minPhi,
            out double XMax, out double XMin, out double YMax, out double YMin)
        { 
            double phi = 0.0, r = 0.0, XC = 0.0, YC = 0.0;

            minPhi = 0.0;
            minR = double.MaxValue;
            XMax = double.MinValue;
            YMax = double.MinValue;
            XMin = double.MaxValue;
            YMin = double.MaxValue;

            while (phi <= 8.0 * Math.PI)
            {
                r = 1.0 / u1(phi);

                if (r < minR)
                {
                    minR = r;
                    minPhi = phi;
                }

                XC = X(r, phi);

                if (XC > XMax)
                    XMax = XC;

                YC = Y(r, phi);

                if (YC > YMax)
                    YMax = YC;

                if (XC < XMin)
                    XMin = XC;

                if (YC < YMin)
                    YMin = YC;

                phi += 0.001;
            }
        }

        private void Panel1_Paint(object sender, PaintEventArgs e)
        {
            panel1.Size = ClientSize;
            int width = ClientSize.Width;
            int height = ClientSize.Height;
            int deltaX = width / 6;
            int deltaY = height / 6;
            int minX = deltaX;
            int maxX = 5 * deltaX;
            int minY = deltaY;
            int maxY = 5 * deltaY;
            double maxPhi, minPhi, maxR, minR;
            double XMax, XMin, YMax, YMin;
            double UMax, UMin, VMax, VMin;

            Maximums(out maxR, out maxPhi, out XMax, out XMin, out YMax, out YMin);
            Minimums(out minR, out minPhi, out UMax, out UMin, out VMax, out VMin);
            
            double slopeX = (maxX - minX) / (XMax - XMin);
            double slopeY = (minY - maxY) / (YMax - YMin);
            double interX = minX - slopeX * XMin;
            double interY = maxY - slopeY * YMin;
            double chi = 0.0, eta = 0.0, phi = 0.0, r = 0.0, x, y;
            Graphics g = e.Graphics;
            Pen bp = new Pen(Color.Black);
            SolidBrush rb = new SolidBrush(Color.Red);

            g.Clip = new Region(new Rectangle(minX, minY, maxX - minX + 1, maxY - minY + 1));

            for (int i = 0; i < 5; i++)
                g.DrawLine(bp, (i + 1) * deltaX, minY, (i + 1) * deltaX, maxY);

            for (int i = 0; i < 5; i++)
                g.DrawLine(bp, minX, (i + 1) * deltaY, maxX, (i + 1) * deltaY);

            while (phi <= 8.0 * Math.PI)
            {
                r = 1.0 / u1(phi);
                x = X(r, phi);
                y = Y(r, phi);
                chi = slopeX * x + interX;
                eta = slopeY * y + interY;
                g.FillEllipse(rb, (float)chi, (float)eta, 1, 1);
                phi += 0.001;
            }
        }
    }
}