Special Relativity

SpecialRelativity.cs

<Query Kind="Program">
  <Reference>&lt;RuntimeDirectory&gt;\System.Numerics.dll</Reference>
  <Namespace>System.Drawing</Namespace>
  <Namespace>System.Drawing.Drawing2D</Namespace>
  <Namespace>System.Drawing.Text</Namespace>
  <Namespace>System.Drawing.Imaging</Namespace>
  <Namespace>System.Runtime.InteropServices</Namespace>
</Query>

const int FrameRate = 30;

void Main()
{
    ConventionalClocks();
    //LightClocks();
    //FlashingRing();
    //Train();
    //Barn();
}

private void Train()
{
    const double LightningDuration = 1.0;
    const double LeadInTime = 2.0;
    const double LeadOutTime = 3.0;

    const double TrainCarLength = 0.015;
    const double TrainCarHeight = 0.02;
    const double TrainCarSpacing = 0.00375;
    const double TrainCars = 10;
    const double TrainVelocity = 0.3;

    const int Frames = 180;
    const double Duration = LeadInTime + LightningDuration + LeadOutTime;
    const double PerFrameTime = Duration / Frames;

    var universe = new Universe(spatialDimensions: 2);

    var S1 = universe.NativeReferenceFrame;

    var S2 = new WorldLine(
        universe,
        new Event(universe.NativeReferenceFrame, 0, 2, 0),
        new Vector(TrainVelocity, 0));

    var scene = new Scene();

    var trainCarWithSpacing = TrainCarLength + TrainCarSpacing;
    var trainLength = trainCarWithSpacing * TrainCars - TrainCarSpacing;
    var halfTrain = trainLength / 2.0;
	var trainTop = TrainCarHeight / 2;

    var renderTrain = new Action<double>(t =>
    {
        for (var c = 0; c < (int)TrainCars; c++)
        {
            var carColor = ColorUtils.FromHsb(c / TrainCars, b: 0.4);
            var x1 = (TrainVelocity * t) + (trainCarWithSpacing * c) - halfTrain;
            var x2 = x1 + TrainCarLength;
            var y1 = trainTop;
            var y2 = -y1;
            var x1y1 = new Event(S1, t, x1, y1);
            var x1y2 = new Event(S1, t, x1, y2);
            var x2y1 = new Event(S1, t, x2, y1);
            var x2y2 = new Event(S1, t, x2, y2);
            scene.Add("", new Segment(x1y1, x1y2), carColor);
            scene.Add("", new Segment(x1y1, x2y1), carColor);
            scene.Add("", new Segment(x2y2, x1y2), carColor);
            scene.Add("", new Segment(x2y2, x2y1), carColor);
        }
    });

    renderTrain(0);
	scene.Add("Left", new Segment(new Event(S1, -LightningDuration, -halfTrain, -(LightningDuration * 0.9) - trainTop), new Event(S1, 0, -halfTrain, -trainTop)), Color.SkyBlue);
	scene.Add("Right", new Segment(new Event(S1, -LightningDuration, halfTrain, -(LightningDuration * 0.9) - trainTop), new Event(S1, 0, halfTrain, -trainTop)), Color.SkyBlue);
	for (var i = 0; ; i++)
    {
		var t = i / (double)FrameRate;
		if (t >= LeadInTime + LightningDuration && t >= LeadOutTime)
		{
			break;
		}

        if (t > 0 && t < LeadInTime + LightningDuration)
        {
            renderTrain(-t);
        }

        if (t < LeadOutTime)
        {
            renderTrain(t);
        }
    }

    RenderPerspectivesToVideo("trains", scene, new[]
    {
        new Scene.Perspective(nameof(S1), S1),
        new Scene.Perspective(nameof(S2), S2),
    }, perFrameTime: PerFrameTime).Dump();
}

private void Barn()
{
    const double PoleLength = 0.00501;
    const double BarnWidth = 0.005;
    const double ExtraContractionFactor = 0.5;
    const double ObservedPoleLength = BarnWidth * ExtraContractionFactor;
    double PoleVelocity = Math.Sqrt(PoleLength * PoleLength - ObservedPoleLength * ObservedPoleLength) / PoleLength;
    double FullLength = ObservedPoleLength + BarnWidth + ObservedPoleLength;
    double Duration = FullLength / PoleVelocity;
    const int Frames = 90;
    double PerFrameTime = Duration / Frames;
    double DoorCloseTime = ObservedPoleLength / PoleVelocity;
    double DoorOpenTime = BarnWidth / PoleVelocity;

    var universe = new Universe(spatialDimensions: 2);

    var S1 = universe.NativeReferenceFrame;

    var S2 = new WorldLine(
        universe,
        new Event(universe.NativeReferenceFrame, 0, 2, 0),
        new Vector(PoleVelocity, 0));

    var scene = new Scene();
    for (var i = 0; i < Frames; i++)
    {
        var t = i * PerFrameTime;

        scene.Add("", new Segment(new Event(S1, t, -BarnWidth / 2, -BarnWidth / 2), new Event(S1, t, BarnWidth / 2, -BarnWidth / 2)), Color.Red);
        scene.Add("", new Segment(new Event(S1, t, BarnWidth / 2, BarnWidth / 2), new Event(S1, t, -BarnWidth / 2, BarnWidth / 2)), Color.Red);
        scene.Add("", new Segment(new Event(S1, t, -BarnWidth / 2, -BarnWidth / 2), new Event(S1, t, 0, -2 * BarnWidth / 3)), Color.Red);
        scene.Add("", new Segment(new Event(S1, t, 0, -2 * BarnWidth / 3), new Event(S1, t, BarnWidth / 2, -BarnWidth / 2)), Color.Red);

        if (t > DoorCloseTime && t < DoorOpenTime)
        {
            scene.Add("", new Segment(new Event(S1, t, BarnWidth / 2, BarnWidth / 2), new Event(S1, t, BarnWidth / 2, -BarnWidth / 2)), Color.Red);
            scene.Add("", new Segment(new Event(S1, t, -BarnWidth / 2, -BarnWidth / 2), new Event(S1, t, -BarnWidth / 2, BarnWidth / 2)), Color.Red);
        }

        var offset = -BarnWidth / 2 - ObservedPoleLength + PoleVelocity * t;
        scene.Add("", new Segment(new Event(S1, t, offset, 0), new Event(S1, t, offset + ObservedPoleLength, 0)), Color.Blue);
    }

    RenderPerspectivesToVideo("barn", scene, new[]
    {
        new Scene.Perspective(nameof(S1), S1),
        new Scene.Perspective(nameof(S2), S2),
    }, perFrameTime: PerFrameTime).Dump();
}

private void LightClocks()
{
    var universe = new Universe(spatialDimensions: 2);

    var S1 = universe.NativeReferenceFrame;

    var S2 = new WorldLine(
        universe,
        new Event(universe.NativeReferenceFrame, 0, 2, 0),
        new Vector(-0.9, 0));

    var scene = new Scene
    {
        { nameof(S1), S1, Color.Orange },
        { nameof(S2), S2, Color.OrangeRed },
    };

    var renderLightClock = new Action<string, WorldLine>((name, o) =>
    {
        for (var i = 0; i < 200; i++)
        {
            var t = i / (double)FrameRate;
            var y = (t - 0.5 - Math.Floor(t)) * Math.Pow(-1, Math.Floor(t));
			scene.Add("", new Segment(new Event(o, t, -0.5, -0.5), new Event(o, t, 0.5, -0.5)), Color.Black);
			scene.Add("", new Segment(new Event(o, t, 0.5, 0.5), new Event(o, t, 0.5, -0.5)), Color.Black);
			scene.Add("", new Segment(new Event(o, t, 0.5, 0.5), new Event(o, t, -0.5, 0.5)), Color.Black);
			scene.Add("", new Segment(new Event(o, t, -0.5, -0.5), new Event(o, t, -0.5, 0.5)), Color.Black);
			scene.Add($"{name} {t:F2}", new Event(o, t, 0, y), ColorUtils.FromHsb(i / 200.0, b: 0.4));
        }
    });
    renderLightClock(nameof(S1), S1);
    renderLightClock(nameof(S2), S2);

    RenderPerspectivesToVideo("light-clocks", scene, new[]
    {
        new Scene.Perspective(nameof(S1), S1),
        new Scene.Perspective(nameof(S2), S2),
    }, fadeFrames: 0).Dump();
}

private void ConventionalClocks()
{
    const double ClockRadius = 1.0;
    const double RevolutionTime = 10.0;
    const double Revolutions = 6.0;

    var universe = new Universe(spatialDimensions: 2);

    var S1 = universe.NativeReferenceFrame;

    var S2 = new WorldLine(
        universe,
        new Event(universe.NativeReferenceFrame, 0, 2, 0),
        new Vector(-0.9, 0));

    var scene = new Scene
    {
        { nameof(S1), S1, Color.Orange },
        { nameof(S2), S2, Color.OrangeRed },
    };

    var renderClock = new Action<string, WorldLine>((name, o) =>
    {
        for (var i = 0; i < FrameRate * RevolutionTime * Revolutions; i++)
        {
            var t = i / (double)FrameRate;
            var x = ClockRadius * Math.Sin(t / RevolutionTime);
            var y = ClockRadius * Math.Cos(t / RevolutionTime);
            scene.Add("", new Segment(new Event(o, t, 0, 0), new Event(o, t, x, y)), Color.Black);
            scene.Add($"{name} {t:F2}", new Event(o, t, x, y), ColorUtils.FromHsb(i / 200.0, b: 0.4));
        }
    });
    renderClock(nameof(S1), S1);
    renderClock(nameof(S2), S2);

    RenderPerspectivesToVideo("conventional-clocks", scene, new[]
    {
        new Scene.Perspective(nameof(S1), S1)
        {
            Focus = MathUtils.LerpXY(new PointF(-2f, -2f), new PointF(0, 0), new PointF(2f, 2f), new PointF(1, 1)),
        },
        new Scene.Perspective(nameof(S2), S2)
        {
            Focus = MathUtils.LerpXY(new PointF(-2f, -2f), new PointF(0, 0), new PointF(2f, 2f), new PointF(1, 1)),
        },
    }, fadeFrames: 0).Dump();
}

private void FlashingRing()
{
    var universe = new Universe(spatialDimensions: 2);

    var S1 = universe.NativeReferenceFrame;

    var S2 = new WorldLine(
        universe,
        new Event(universe.NativeReferenceFrame, 0, 0.15, 0),
        new Vector(-0.99, 0));

    var S3 = new WorldLine(
        universe,
        new Event(universe.NativeReferenceFrame, 0, 2, 0),
        new Vector(-0.99, 0));

    var scene = new Scene
    {
        { nameof(S1), S1, Color.Orange },
        { nameof(S2), S2, Color.OrangeRed },
        { nameof(S3), S3, Color.DarkOrange },
    };

    for (var i = 0; i < 100; i++)
    {
        var r = i / 100.0;
        scene.Add(i.ToString("00"), new Event(S1, 0, Math.Sin(r * MathUtils.TAU), Math.Cos(r * MathUtils.TAU)), ColorUtils.FromHsb(r, b: 0.4));
    }

    RenderPerspectivesToVideo("flashing-ring", scene, new[]
    {
        new Scene.Perspective(nameof(S1), S1)
        {
            Focus = MathUtils.LerpXY(new PointF(-2.3f, -2.3f), new PointF(0, 0), new PointF(2.3f, 2.3f), new PointF(1, 1)),
        },
        new Scene.Perspective(nameof(S2), S2),
        new Scene.Perspective(nameof(S3), S3),
    }).Dump();
}

private object RenderPerspectivesToVideo(string title, Scene scene, Scene.Perspective[] perspectives, double perFrameTime = 1D / FrameRate, int width = 300, int margin = 20, int fadeFrames = 6, int leadFrames = 30)
{
    var correctedFrames = Array.ConvertAll(perspectives, p => scene.Render(p, width, perFrameTime, fadeFrames, abberation: false).ToList());
    var flightTimeFrames = Array.ConvertAll(perspectives, p => scene.Render(p, width, perFrameTime, fadeFrames, abberation: true).ToList());
    var maxFrames = Math.Max(correctedFrames.Max(f => f.Count), flightTimeFrames.Max(f => f.Count));
    var outputPath = Path.Combine(Environment.GetFolderPath(Environment.SpecialFolder.Desktop), $"{title}.gif");
    return Enumerable
        .Range(0, maxFrames + leadFrames)
        .Select(i => ImageUtils.Render(margin + (width + margin) * perspectives.Length, margin + (margin + width) * 2, g =>
        {
            g.Clear(Color.White);
            for (var j = 0; j < perspectives.Length; j++)
            {
                var topLeft = new Point(margin + j * (width + margin), margin);
                var middleLeft = new Point(margin + j * (width + margin), margin + (width + margin));
                if (i >= leadFrames)
                {
                    var f = i - leadFrames;
                    var cf = correctedFrames[j];
                    if (f < cf.Count)
                    {
                        g.DrawImageUnscaled(cf[f], topLeft);
                        cf[f].Dispose();
                    }

                    var ff = flightTimeFrames[j];
                    if (f < ff.Count)
                    {
                        g.DrawImageUnscaled(ff[f], middleLeft);
                        ff[f].Dispose();
                    }
                }

                g.DrawRectangle(Pens.Black, new Rectangle(topLeft, new Size(width, width)));
                g.DrawRectangle(Pens.Black, new Rectangle(middleLeft, new Size(width, width)));
            }
        }))
        .TrackProgress(maxFrames + leadFrames, outputPath)
        .ToAnimation(outputPath, FrameRate);
}

class Scene : IEnumerable<Scene.ITagged>
{
    private readonly List<ITagged> elements = new List<ITagged>();

    public void Add<T>(string name, T value, Color color)
    {
        this.elements.Add(new Tagged<T>(name, value, color));
    }

    public IEnumerator<ITagged> GetEnumerator() => this.elements.GetEnumerator();

    IEnumerator IEnumerable.GetEnumerator() => this.GetEnumerator();

    public IEnumerable<Bitmap> Render(Perspective perspective, int width, double perFrameTime = 1D / FrameRate, int fadeFrames = 6, bool abberation = false)
    {
        const float eventDrawRadius = 5.0f;
		var font = SystemFonts.DefaultFont;
		var withAbberation = new Func<Event, Event>(e => abberation ? new Event(e.ReferenceFrame, e[0] + Math.Sqrt(e[1] * e[1] + e[2] * e[2]) / e.ReferenceFrame.Universe.C, e[1], e[2]) : e);

		var worldLines = this.elements.OfType<ITagged<WorldLine>>();
        var allSegments = this.elements.OfType<ITagged<Segment>>();
        var allEvents = worldLines.Select(w => new Tagged<Event>(w.Name, w.Value.Origin, w.Color)).Concat(this.elements.OfType<ITagged<Event>>());
        var mappedSegments = (from s in allSegments
                              let convertedA = withAbberation(s.Value.EventA.ToReferenceFrame(perspective.ReferenceFrame))
                              let convertedB = withAbberation(s.Value.EventB.ToReferenceFrame(perspective.ReferenceFrame))
                              let newSegment = new Segment(convertedA, convertedB)
                              orderby Math.Min(convertedA[0], convertedB[0])
                              select new Tagged<Segment>(s.Name, newSegment, s.Color)).ToList();
        var mappedEvents = (from e in allEvents
                            let converted = withAbberation(e.Value.ToReferenceFrame(perspective.ReferenceFrame))
                            orderby converted[0]
                            select new Tagged<Event>(e.Name, converted, e.Color)).ToList();

        var extents =
            mappedEvents.Select(e => e.Value)
            .Concat(mappedSegments.SelectMany(s => new[] { s.Value.EventA, s.Value.EventB }))
            .Aggregate(new
            {
                min0 = double.NaN,
                max0 = double.NaN,
                min1 = double.NaN,
                max1 = double.NaN,
                min2 = double.NaN,
                max2 = double.NaN,
            },
            (a, b) => new
            {
                min0 = double.IsNaN(a.min0) ? b[0] : Math.Min(a.min0, b[0]),
                max0 = double.IsNaN(a.max0) ? b[0] : Math.Max(a.max0, b[0]),
                min1 = double.IsNaN(a.min1) ? b[1] : Math.Min(a.min1, b[1]),
                max1 = double.IsNaN(a.max1) ? b[1] : Math.Max(a.max1, b[1]),
                min2 = double.IsNaN(a.min2) ? b[2] : Math.Min(a.min2, b[2]),
                max2 = double.IsNaN(a.max2) ? b[2] : Math.Max(a.max2, b[2]),
            });
        var frameTime = new Func<int, double>(f => extents.min0 + f * perFrameTime);
        var frames = Math.Max((int)Math.Ceiling((extents.max0 - extents.min0) / perFrameTime), 1) + fadeFrames;

        var frameEvents = new List<Tagged<Event>>[frames];
        var eventEndIndex = 0;
        for (var f = 0; f < frames; f++)
        {
            var events = frameEvents[f] = new List<Tagged<Event>>();
            var eventStartIndex = eventEndIndex;
            var endTime = frameTime(f + 1);

            while (eventEndIndex < mappedEvents.Count && mappedEvents[eventEndIndex].Value[0] < endTime)
            {
                eventEndIndex++;
            }

            var eventCount = eventEndIndex - eventStartIndex;
            if (eventCount > 0)
            {
                events.AddRange(mappedEvents.GetRange(eventStartIndex, eventCount));
            }
        }

        var frameSegments = new List<Tagged<Segment>>[frames];
        for (var f = 0; f < frames; f++)
        {
            var startTime = frameTime(f);
            var endTime = frameTime(f + 1);
            frameSegments[f] = mappedSegments.Where(s => s.Value.EventB[0] >= startTime && s.Value.EventA[0] < endTime || s.Value.EventA[0] >= startTime && s.Value.EventB[0] < endTime).ToList();
        }

        var focus = perspective.Focus;
        if (focus == null)
        {
            var diffX = extents.max1 - extents.min1;
            var diffY = extents.max2 - extents.min2;
            if (diffX == 0) diffX = 10;
            if (diffY == 0) diffY = 10;
            var left = extents.min1 - diffX / 10;
            var right = extents.max1 + diffX / 10;
            var top = extents.min2 - diffY / 10;
            var bottom = extents.max2 + diffY / 10;
            focus = MathUtils.LerpXY(new PointF((float)left, (float)top), new PointF(0, 0), new PointF((float)right, (float)bottom), new PointF(1, 1));
        }

        for (var f = 0; f < frames; f++)
        {
            var frameStart = frameTime(f);
            var frameEnd = frameTime(f + 1);
            yield return ImageUtils.Render(width, width, g =>
            {
                g.Clear(Color.White);

                g.InterpolationMode = InterpolationMode.NearestNeighbor;
                g.SmoothingMode = SmoothingMode.None;
                g.TextRenderingHint = TextRenderingHint.SingleBitPerPixelGridFit;

                //g.InterpolationMode = InterpolationMode.High;
                //g.SmoothingMode = SmoothingMode.HighQuality;
                //g.TextRenderingHint = TextRenderingHint.AntiAlias;

                var caption = perspective.Name;
                if (abberation)
                {
                    caption += " (abberation)";
                }

                if (true)
                {
                    var precision = (int)Math.Ceiling(Math.Log10(frames) - Math.Log10(extents.max0 - extents.min0));
                    caption += " t = " + frameStart.ToString("F" + precision);
                }

                var size = g.MeasureString(caption, font);
                g.DrawString(caption, font, Brushes.Black, 0, width - size.Height);

                var drawEvent = new Action<ITagged<Event>, double>((frameEvent, fade) =>
                {
                    var c = focus(new PointF((float)frameEvent.Value[1], (float)frameEvent.Value[2]));
                    var x = c.X * width;
                    var y = c.Y * width;
                    using (var brush = new SolidBrush(ColorUtils.Blend(fade, Color.White, frameEvent.Color)))
                    {
                        g.FillEllipse(brush, new RectangleF(x - eventDrawRadius / 2, y - eventDrawRadius / 2, eventDrawRadius, eventDrawRadius));
                        if (!string.IsNullOrEmpty(frameEvent.Name))
                        {
                            g.DrawString(frameEvent.Name, font, brush, x + eventDrawRadius, y + eventDrawRadius);
                        }
                    }
                });

                var drawSegment = new Action<ITagged<Segment>>(frameSegment =>
                {
                    var eventA = frameSegment.Value.EventA;
                    var eventB = frameSegment.Value.EventB;
                    var a0 = eventA[0];
                    var b0 = eventB[0];
                    var a1 = eventA[1];
                    var b1 = eventB[1];
                    var a2 = eventA[2];
                    var b2 = eventB[2];

                    if (a0 < frameStart || a0 >= frameEnd ||
                        b0 < frameStart || b0 >= frameEnd)
                    {
                        var a0_ = MathUtils.Clamp(a0, frameStart, frameEnd);
                        var b0_ = MathUtils.Clamp(b0, frameStart, frameEnd);

                        var lerp1 = MathUtils.Lerp(a0, a1, b0, b1);
                        var lerp2 = MathUtils.Lerp(a0, a2, b0, b2);

                        if (a0_ != a0)
						{
							a0 = a0_;
							a1 = lerp1(a0_);
							a2 = lerp2(a0_);
						}

                        if (b0_ != b0)
						{
							b0 = b0_;
							b1 = lerp1(b0_);
							b2 = lerp2(b0_);
						}
                    }

                    var cA = focus(new PointF((float)a1, (float)a2));
                    var cB = focus(new PointF((float)b1, (float)b2));
                    var xA = cA.X * width;
                    var yA = cA.Y * width;
                    var xB = cB.X * width;
                    var yB = cB.Y * width;
                    using (var pen = new Pen(frameSegment.Color))
                    {
                        g.DrawLine(pen, xA, yA, xB, yB);
                        if (!string.IsNullOrEmpty(frameSegment.Name))
                        {
                            using (var brush = new SolidBrush(frameSegment.Color))
                            {
                                //g.DrawString(frameSegment.Name, font, brush, x, y);
                            }
                        }
                    }
                });

                for (var fade = fadeFrames; fade >= 1; fade--)
                {
                    if (f - fade < 0)
                    {
                        continue;
                    }

                    foreach (var frameEvent in frameEvents[f - fade])
                    {
                        drawEvent(frameEvent, (double)fade / fadeFrames);
                    }
                }

                foreach (var frameEvent in frameEvents[f])
                {
                   drawEvent(frameEvent, 0);
                }
                
                foreach (var frameSegment in frameSegments[f])
                {
                    drawSegment(frameSegment);
                }
            });
        }
    }

    public interface ITagged
    {
        string Name { get; }
        Color Color { get; }
    }

    public interface ITagged<out T> : ITagged
    {
        T Value { get; }
    }

    public class Perspective
    {
        public Perspective(string name, WorldLine referenceFrame)
        {
            this.Name = name;
            this.ReferenceFrame = referenceFrame;
        }

        public string Name { get; set; }
        public WorldLine ReferenceFrame { get; set; }
        public Func<PointF, PointF> Focus { get; set; }
    }

    private class Tagged<T> : ITagged<T>
    {
        public Tagged(string name, T value, Color color)
        {
            this.Name = name;
            this.Value = value;
            this.Color = color;
        }

        public string Name { get; set; }
        public T Value { get; set; }
        public Color Color { get; set; }
    }
}

class Vector
{
    private readonly double[] coordinates;

    public Vector(params double[] coordinates)
    {
        if (object.ReferenceEquals(coordinates, null))
        {
            throw new ArgumentNullException(nameof(coordinates));
        }
        else if (coordinates.Length < 1)
        {
            throw new ArgumentOutOfRangeException(nameof(coordinates));
        }

        this.coordinates = coordinates.ToArray();
    }

    public int Dimensions => this.coordinates.Length;

    public double this[int index] => this.coordinates[index];

    public double Magnitude => Math.Sqrt(Vector.DotProduct(this, this));

    public static double DotProduct(Vector a, Vector b)
    {
        if (object.ReferenceEquals(a, null))
        {
            throw new ArgumentNullException(nameof(a));
        }
        else if (object.ReferenceEquals(b, null))
        {
            throw new ArgumentNullException(nameof(b));
        }
        else if (a.Dimensions != b.Dimensions)
        {
            throw new ArgumentOutOfRangeException(nameof(b));
        }

        return Matrix.MultiplyRowColumn(a.coordinates, b.coordinates);
    }

    public static Vector operator -(Vector operand)
    {
        if (object.ReferenceEquals(operand, null))
        {
            throw new ArgumentNullException(nameof(operand));
        }

        return new Vector(operand.coordinates.Select(c => -c).ToArray());
    }

    public static Vector operator *(double left, Vector right)
    {
        if (object.ReferenceEquals(right, null))
        {
            throw new ArgumentNullException(nameof(right));
        }

        return new Vector(right.coordinates.Select(c => left * c).ToArray());
    }

    public override string ToString() => $"[{string.Join(", ", this.coordinates)}]";
}

class Event
{
    private readonly double[] coordinates;

    public Event(WorldLine referenceFrame, params double[] coordinates)
    {
        if (object.ReferenceEquals(referenceFrame, null))
        {
            throw new ArgumentNullException(nameof(referenceFrame));
        }
        else if (object.ReferenceEquals(coordinates, null))
        {
            throw new ArgumentNullException(nameof(coordinates));
        }
        else if (coordinates.Length != referenceFrame.Universe.Dimensions)
        {
            throw new ArgumentOutOfRangeException(nameof(coordinates));
        }

        this.ReferenceFrame = referenceFrame;        
        this.coordinates = coordinates.ToArray();
    }

    public WorldLine ReferenceFrame { get; }

    public double this[int index] => this.coordinates[index];

    public static double IntervalSq(Event a, Event b)
    {
        if (object.ReferenceEquals(a, null))
        {
            throw new ArgumentNullException(nameof(a));
        }
        else if (object.ReferenceEquals(b, null))
        {
            throw new ArgumentNullException(nameof(b));
        }
        else if (!object.ReferenceEquals(a.ReferenceFrame, b.ReferenceFrame))
        {
            throw new ArgumentException($"Reference frame of '{nameof(b)}' must match the reference frame of '{nameof(a)}'.", nameof(b));
        }

        var cSq = a.ReferenceFrame.Universe.C;
        return Enumerable
            .Range(0, a.ReferenceFrame.Universe.SpatialDimensions + 1)
            .Select(i => b[i] - a[i])
            .Select((v, i) => (i == 0 ? -1 / cSq : 1) * (v * v))
            .Sum();
    }

    public static double[,] LorentzTransform(Universe universe, Vector relativeVelocity)
    {
        var d = universe.Dimensions;
        var v = relativeVelocity;
        var c = universe.C;
        var cSq = c * c;
        var vSq = Vector.DotProduct(v, v);
        var gamma = 1 / Math.Sqrt(1 - vSq / cSq);

        var transform = new double[d, d];
        transform[0, 0] = gamma;
        for (var a = 1; a < d; a++)
        {
            transform[0, a] = transform[a, 0] = -(v[a - 1] / c) * gamma;
            transform[a, a] = (gamma - 1) * ((v[a - 1] * v[a - 1]) / vSq) + 1;
            for (var b = a + 1; b < d; b++)
            {
                transform[a, b] = transform[b, a] = (gamma - 1) * ((v[a - 1] * v[b - 1]) / vSq);
            }
        }

        return transform;
    }

    public Event ToNativeReferenceFrame()
    {
        if (object.ReferenceEquals(this.ReferenceFrame.Universe.NativeReferenceFrame, this.ReferenceFrame))
        {
            return this;
        }

        var srcFrame = this.ReferenceFrame;
        var universe = srcFrame.Universe;
        var d = universe.Dimensions;

        var transform = Event.LorentzTransform(universe, srcFrame.RelativeVelocity);
        var boostedOrigin = Matrix.MultiplyColumn(transform, srcFrame.Origin.coordinates);

        var boosted = new double[d];
        for (var i = 0; i < d; i++)
        {
            boosted[i] = this.coordinates[i] + boostedOrigin[i];
        }

        for (var a = 1; a < d; a++)
        {
            transform[a, 0] = transform[0, a] = -transform[a, 0];
        }

        var native = Matrix.MultiplyColumn(transform, boosted);

        return new Event(universe.NativeReferenceFrame, native);
    }

    public Event ToReferenceFrame(WorldLine referenceFrame)
    {
        if (object.ReferenceEquals(referenceFrame, null))
        {
            throw new ArgumentNullException(nameof(referenceFrame));
        }
        else if (object.ReferenceEquals(referenceFrame, this.ReferenceFrame))
        {
            return this;
        }
        else if (!object.ReferenceEquals(this.ReferenceFrame.Universe, referenceFrame.Universe))
        {
            throw new ArgumentException($"Universe of '{nameof(referenceFrame)}' must match the universe of the source event's reference frame.", nameof(referenceFrame));
        }

        var universe = this.ReferenceFrame.Universe;
        var d = universe.Dimensions;
        var nativeEvent = this.ToNativeReferenceFrame();

        if (referenceFrame == nativeEvent.ReferenceFrame)
        {
            return nativeEvent;
        }

        var transform = Event.LorentzTransform(universe, referenceFrame.RelativeVelocity);
        var boosted = Matrix.MultiplyColumn(transform, nativeEvent.coordinates);
        var boostedOrigin = Matrix.MultiplyColumn(transform, referenceFrame.Origin.coordinates);

        var translated = new double[d];
        for (var i = 0; i < d; i++)
        {
            translated[i] = boosted[i] - boostedOrigin[i];
        }

        return new Event(referenceFrame, translated);
    }

    public override string ToString() => $"({string.Join(", ", this.coordinates)})";
}

class Segment
{
    public Segment(Event eventA, Event eventB)
    {
        if (object.ReferenceEquals(eventA, null))
        {
            throw new ArgumentNullException(nameof(eventA));
        }
        else if (object.ReferenceEquals(eventB, null))
        {
            throw new ArgumentNullException(nameof(eventB));
        }

        this.EventA = eventA;
        this.EventB = eventB.ToReferenceFrame(eventA.ReferenceFrame);
    }

    public Event EventA { get; }

    public Event EventB { get; }
}

class Universe
{
    public Universe(double c = 1, int spatialDimensions = 3)
    {
        if (double.IsNaN(c) || double.IsInfinity(c))
        {
            throw new ArgumentOutOfRangeException(nameof(c));
        }
        else if (spatialDimensions < 1)
        {
            throw new ArgumentOutOfRangeException(nameof(spatialDimensions));
        }

        this.C = c;
        this.SpatialDimensions = spatialDimensions;
        this.NativeReferenceFrame = new WorldLine(this);
    }

    public double C { get; }
    public WorldLine NativeReferenceFrame { get; }
    public int SpatialDimensions { get; }
    public int Dimensions => this.SpatialDimensions + 1;
}

class WorldLine
{
    internal WorldLine(Universe universe)
    {
        this.Universe = universe;
        this.Origin = new Event(this, new double[universe.Dimensions]);
        this.RelativeVelocity = new Vector(new double[universe.SpatialDimensions]);
    }

    public WorldLine(Universe universe, Event origin, Vector relativeVelocity)
    {
        if (object.ReferenceEquals(universe, null))
        {
            throw new ArgumentNullException(nameof(universe));
        }
        else if (object.ReferenceEquals(origin, null))
        {
            throw new ArgumentNullException(nameof(origin));
        }
        else if (!object.ReferenceEquals(universe.NativeReferenceFrame, origin.ReferenceFrame))
        {
            throw new ArgumentException($"Reference frame of '{nameof(origin)}' must match the native reference frame of '{nameof(universe)}'.", nameof(origin));
        }
        else if (object.ReferenceEquals(relativeVelocity, null))
        {
            throw new ArgumentNullException(nameof(relativeVelocity));
        }
        else if (relativeVelocity.Dimensions != universe.SpatialDimensions)
        {
            throw new ArgumentOutOfRangeException(nameof(relativeVelocity));
        }

        this.Universe = universe;
        this.Origin = origin;
        this.RelativeVelocity = relativeVelocity;
    }

    public Universe Universe { get; }
    public Event Origin { get; }
    public Vector RelativeVelocity { get; }

    public override string ToString() => $"{this.Origin} -> {this.RelativeVelocity}";
}

public static class Matrix
{
    public static double[,] Multiply(double[,] a, double[,] b)
    {
        if (object.ReferenceEquals(a, null))
        {
            throw new ArgumentNullException(nameof(a));
        }
        else if (object.ReferenceEquals(b, null))
        {
            throw new ArgumentNullException(nameof(b));
        }
        else if (a.GetLength(1) != b.GetLength(0))
        {
            throw new ArgumentOutOfRangeException(nameof(b));
        }

        var n = a.GetLength(0);
        var m = a.GetLength(1);
        var p = b.GetLength(1);

        var result = new double[n, p];
        for (var i = 0; i < n; i++)
        {
            for (var j = 0; j < p; j++)
            {
                for (var k = 0; k < m; k++)
                {
                    result[i, j] += a[i, k] * b[k, j];
                }
            }
        }

        return result;
    }

    public static double[] MultiplyColumn(double[,] a, double[] b)
    {
        if (object.ReferenceEquals(a, null))
        {
            throw new ArgumentNullException(nameof(a));
        }
        else if (object.ReferenceEquals(b, null))
        {
            throw new ArgumentNullException(nameof(b));
        }
        else if (a.GetLength(1) != b.Length)
        {
            throw new ArgumentOutOfRangeException(nameof(b));
        }

        var n = a.GetLength(0);
        var m = a.GetLength(1);

        var result = new double[n];
        for (var i = 0; i < n; i++)
        {
            for (var k = 0; k < m; k++)
            {
                result[i] += a[i, k] * b[k];
            }
        }

        return result;
    }

    public static double[] MultiplyRow(double[] a, double[,] b)
    {
        if (object.ReferenceEquals(a, null))
        {
            throw new ArgumentNullException(nameof(a));
        }
        else
        if (object.ReferenceEquals(b, null))
        {
            throw new ArgumentNullException(nameof(b));
        }
        else if (a.Length != b.GetLength(0))
        {
            throw new ArgumentOutOfRangeException(nameof(b));
        }

        var m = a.Length;
        var p = b.GetLength(1);

        var result = new double[p];
        for (var j = 0; j < p; j++)
        {
            for (var k = 0; k < m; k++)
            {
                result[j] += a[k] * b[k, j];
            }
        }

        return result;
    }

    public static double MultiplyRowColumn(double[] a, double[] b)
    {
        if (object.ReferenceEquals(a, null))
        {
            throw new ArgumentNullException(nameof(a));
        }
        else if (object.ReferenceEquals(b, null))
        {
            throw new ArgumentNullException(nameof(b));
        }
        else if (a.Length != b.Length)
        {
            throw new ArgumentOutOfRangeException(nameof(b));
        }

        var m = a.Length;

        double result = default(double);
        for (var k = 0; k < m; k++)
        {
            result += a[k] * b[k];
        }

        return result;
    }

    public static double[,] MultiplyColumnRow(double[] a, double[] b)
    {
        if (object.ReferenceEquals(a, null))
        {
            throw new ArgumentNullException(nameof(a));
        }
        else if (object.ReferenceEquals(b, null))
        {
            throw new ArgumentNullException(nameof(b));
        }

        var n = a.Length;
        var p = b.Length;

        var result = new double[n, p];
        for (var i = 0; i < n; i++)
        {
            for (var j = 0; j < p; j++)
            {
                result[i, j] += a[i] * b[j];
            }
        }

        return result;
	}
}