thedeemon | В искривленном пространстве. Часть 2

Ну и вот. Написал я тогда модуль символьных вычислений, который умеет арифметику, умеет частные производные, в том числе тригонометрии и сложных ф-й, умеет немного упрощать (выносить за скобки, сокращать дроби, знает, что sin^2 + cos^2 = 1) и умеет генерить код. Этот модуль используется в компайл-тайме, и из одного лишь уравнения поверхности выводит все необходимые формулы и генерит код для каждой фигуры, который тут же компилятором оптимизируется. Посмотрим на примере все той же сферы.
Задаем уравнение поверхности, отображение 2D координат (u,v) в 3D (x,y,z).

Expr[] sphereEq() {
    auto R = new Var("R");
    return [mul(R, mul(new Cos("u"), new Cos("v"))),
            mul(R, new Sin("v")),
            mul(R, mul(new Sin("u"), new Cos("v")))  ];
}

Ф-я возвращает массив из трех выражений - для x, для y и для z. X у меня идет вправо, Y вверх и Z вдаль. Координата u это долгота, координата v - широта, считается от экватора. Имея эти выражения как три компоненты ф-ии
X(u,v) = [R * cos(u) * cos(v), R * sin(v), R * cos(v) * sin(u)]
первым делом находим базисные вектора - производные X(u,v) по u и по v:

    auto Xu = X.diff("u"); 
    auto Xv = X.diff("v");

Получаем вектора
Xu: [[-1 * R * cos(v) * sin(u)], 0, [R * cos(u) * cos(v)]]
Xv: [[-1 * R * cos(u) * sin(v)], [R * cos(v)], [-1 * R * sin(u) * sin(v)]]

теперь нам нужны их попарные скалярные произведения, из которых составляется матрица метрического тензора, и тут же находим к ней обратную матрицу. Все символическими вычислениями, без конкретных чисел.

Expr[] diff(Expr[] vec, string dx) { return vec.amap!(e => e.diff(dx)); }
Expr dot(Expr[] a, Expr[] b) { return iota(a.length).amap!(i => mul(a[i], b[i])).add; }


    Expr[2][2] g, ginv; // metric tensor and its inverse
    g[0][0] = dot(Xu, Xu);
    g[0][1] = g[1][0] = dot(Xu, Xv);
    g[1][1] = dot(Xv, Xv);
    auto det_g = add(mul(g[0][0], g[1][1]), neg(mul(g[0][1], g[1][0])));
    ginv[0][0] = div(g[1][1], det_g);
    ginv[0][1] = ginv[1][0] = div(neg(g[0][1]), det_g);
    ginv[1][1] = div(g[0][0], det_g);

Для сферы получим

g[0][0] = R*R*cos(v)*cos(v)  g[0][1] = 0
g[1][0] = 0                  g[1][1] = R*R

Тут видно, что один градус широты имеет одну и ту же длину везде, а вот один градус долготы меняет длину пропорционально косинусу широты.
Имея выражения для метрики и умея их продифференцировать символьно, получаем символы Кристоффеля:

    Expr[2][2][2] C; // Christoffel symbols
    string[] dx = ["u","v"];
    auto two = new Const("2");
    foreach(k; 0..2)
        foreach(i; 0..2)
            foreach(j; 0..2) {
                //Ckij = 1/2 * ginv_kl * (d/dx_j g_li  + d/dx_i g_lj - d/dx_l g_ij )
                Expr[2] tmp;
                foreach(l; 0..2) {
                    auto J = g[l][i].diff(dx[j]);
                    auto I = g[l][j].diff(dx[i]);
                    auto L = neg(g[i][j].diff(dx[l]));
                    tmp[l] = div( mul(ginv[k][l], add(J, I, L)), two);
                }
                C[k][i][j] = add(tmp[0], tmp[1]);
            }

Для сферы это
C[0][0][1] = C[0][1][0] = -sin(v) / cos(v)
C[1][0][0] = cos(v) * sin(v)
остальные все 0.

Дальше в рантайме нам надо уметь решать уравнение геодезической
d2 x_j + Г_jik * dx_k * dx_i = 0
где d и d2 - первая и вторая производные по абстрактному параметру λ, по i и k идет суммирование, а по j получаются два независимых уравнения - по одному на координату. Т.е. для координат u и v получаем два уравнения
ddu + Г_000 * du * du + 2*Г_001 * du * dv + Г_011 * dv * dv = 0
ddv + Г_100 * du * du + 2*Г_101 * du * dv + Г_111 * dv * dv = 0

Решалку диффуров первой степени я делал пару лет назад, так что достаточно было перевести на D:

// f'(t,y) = drv(y)
T evolveRK(alias drv, T)(T state, double dt) {
    auto a = drv(state);
    auto b = drv(state + a * (dt/2));
    auto c = drv(state + b * (dt/2));
    auto d = drv(state + c*dt);
    return state + a*(dt/6) + b*(dt/3) + c*(dt/3) + d*(dt/6);
}

Тут Т - произвольный тип, для которого есть операции сложения и умножения на число. Например, вектор из даблов.
Тогда если у нас state это текущие значения координат и скоростей (u, du, v, dv), то их производные по λ это будет набор state' = (du, ddu, dv, ddv). Если у нас есть ф-я, вычисляющая state' из state, то ее можно передать в решалку диффуров и получить ответ - последовательность значений state, траекторию в координатах u и v, той же широте и долготе. А получить (du, ddu, dv, ddv) из (u, du, v, dv) понятно как - du и dv просто копируем как есть, а ddu и ddv вычисляем из нашего уравнения выше:
ddu = - Г_000 * du * du - 2*Г_001 * du * dv - Г_011 * dv * dv
ddv = - Г_100 * du * du - 2*Г_101 * du * dv - Г_111 * dv * dv

Один вызов evolveRK дает нам один шаг по геодезической за "время" dt. Длину этого шага можно вычислить как длину вектора (du*dt, dv*dt), используя значения метрического тензора в точке (u,v). Остается лишь делать такие шажки в цикле, пока пройденное расстояние не окажется искомым. Для каждого вертикального столбца пикселей мы один раз проходим по соответствующей геодезической: сперва до места, куда попадает луч через самый нижний пиксель, потом идем дальше до расстояния, где упал второй луч, потом до точки падения третьего и т.д.

Что получилось. Лучше всего это видеть вживую, интерактивно, или хотя бы в движении на видео.

Watch on YouTube

Бинарники под Линукс и Винду можно взять тут. Исходники и краткое описание выложены там же. Виндовая версия примерно вдвое быстрее линуксовой. Даже будучи запущена в линуксе под вайном.

Про результаты и их осмысление я еще буду отдельно писать, а пока кратко.
Кривизна бывает нулевая, положительная и отрицательная. Нулевая - плоское пространство, как мы привыкли. При приближении к объекту он визуально растет, соизмеримо с расстоянием. Сумма углов в треугольнике 180 градусов.
У сферы кривизна положительна, сумма углов больше 180. Лучи из глаз (как завещал Аристотель, просто так проще говорить) сперва расходятся друг от друга, но потом, пройдя 90 градусов сферы, начинают сходиться опять. Поэтому вблизи все выглдядит довольно привычно и обычно, но удаленные объекты ведут себя странно - при приближении могут визуально уменьшаться. И да, если смотреть далеко, то лучи собираются в одном месте, это место визуально бесконечно большое - видно во всех направлениях, а дальше лучи перекрещиваются и лево и право меняются местами.
Хорошо заметно, как радиус сферы влияет на визуальные эффекты: чем больше радиус, тем меньше кривизна, тем ближе картинка к привычной нам.

Сферу можно сплющить, тогда кривизна получается неодинаковой в разных местах. На полюсах такая фигура выглядит более плоской, привычной, а на экваторе наоборот визуальные эффекты кривизны становятся сильнее.

Бублик хорошая фигура. Наружняя его часть похожа на то, что на экваторе эллипсоида, а вот внутренняя часть имеет отрицательную кривизну. В таких местах сумма углов в треугольнике должна быть меньше 180. Лучи из глаз расходятся в стороны и разбегаются все быстрее, нелинейно. Поэтому удаленный объект при приближении к нему может увеличиваться визуально сильнее обычного, часто примерно так же, как более близкий объект. Поэтому когда в ту сторону идешь, вся картинка просто растягивается вширь, будто нарисована на стене. Но могут быть и более хитрые эффекты. Еще на бублике видно, что внутренняя часть меньше внешней, но при этом простраство замкнуто во все стороны, визуально бесконечно, и как такая бесконечная плоскость сочетается с разными длинами ее частей - это надо видеть.

Черная дыра. Про нее еще будет отдельно. Берем метрику Шварцшильда, выкидываем время и одну из пространственных координат, остается плоскость с ЧД в центре. Введя на ней полярные координаты, оказалось несложно найти такую ф-ю высоты Y(r), что выводимая описанным в посте способом метрика совпадает с Шварцшильдовской. Так что все пространственные эффекты искривления пространства там можно пронаблюдать: как оно почти плоское вдали, как удлиняется ближе к ЧД, как свет огибает ЧД и даже наматывает круги вокруг нее. Становится понятно, почему при том, что на горизонте ЧД у нас метрика имеет сингулярность (производная Y(r) становится бесконечной), расстояние до горизонта тем не менее конечно, и не так уж велико.

Ну и wormhole, это надо испытать, непередаваемо. Два бесконечных мира, практически две плоскости, но в одном месте они очень ловко соединяются, можно перейти из одного в другой.