std::array<dual, 4> symmetric_warp_drive(dual t, dual r, dual theta, dual phi)
{
    dual rg = 1;
    dual rk = rg;

    dual a20 = (1 - rg / r);

    dual a0 = sqrt(a20);

    dual a1 = t / theta;

    dual a2 = a20 + a1;

    dual yrr0 = 1 / (1 - (rg / r));
    dual ythetatheta0 = r * r;
    dual yphiphi0 = r * r * sin(theta) * sin(theta);

    dual gamma_0 = pow(r, 4) * sin(theta) * sin(theta) / (1 - rg/r);

    dual littlea = rk * theta * pow(a0, -1);
    dual littleb = rk * theta - sqrt(gamma_0);

    dual Urt = (littlea * pow(a20 + t/theta, 3/2.f) - littleb) / (littlea * a0*a0*a0 - littleb);

    dual yrr = Urt * yrr0;
    dual ythetatheta = Urt * ythetatheta0;
    dual yphiphi = Urt * yphiphi0;

    dual dt = -a2;

    return {dt, yrr, ythetatheta, yphiphi};
}