vec4 lorentzTransform(in vec3 v, in vec4 p)
{
    //https://en.wikipedia.org/wiki/Lorentz_transformation

    const float epsilon = 1.0e-7;

//  float v2 = v.x*v.x + v.y*v.y + v.z*v.z;
    float v2 = dot(v,v);

    float g = 1. / sqrt(1. - v2);
    float gm1 = g - 1.;
    if (v2 < epsilon)
        return p;

    vec4 r;

/*	ORIGINAL CODE
    r.w = g*p.w - v.x*g*p.x - v.y*g*p.y - v.z*g*p.z;
    r.x = -v.x*g*p.w + (1.0 + gm1*v.x*v.x/v2)*p.x + (gm1*v.x*v.y/v2)*p.y + (gm1*v.x*v.z/v2)*p.z;
    r.y = -v.y*g*p.w + (gm1*v.y*v.x/v2)*p.x + (1.0 + gm1*v.y*v.y/v2)*p.y + (gm1*v.y*v.z/v2)*p.z,
    r.z = -v.z*g*p.w + (gm1*v.z*v.x/v2)*p.x + (gm1*v.z*v.y/v2)*p.y + (1.0+gm1*v.z*v.z/v2)*p.z;
*/

/*	STEP 1
    r.w = g * (p.w - (v.x*p.x + v.y*p.y + v.z*p.z));
    r.x = -v.x*g*p.w + p.x + (gm1*v.x*v.x/v2)*p.x + (gm1*v.x*v.y/v2)*p.y + (gm1*v.x*v.z/v2)*p.z;
    r.y = -v.y*g*p.w + p.y + (gm1*v.y*v.x/v2)*p.x + (gm1*v.y*v.y/v2)*p.y + (gm1*v.y*v.z/v2)*p.z,
    r.z = -v.z*g*p.w + p.z + (gm1*v.z*v.x/v2)*p.x + (gm1*v.z*v.y/v2)*p.y + (gm1*v.z*v.z/v2)*p.z;
*/

/*	STEP 2
    vec3 p3 = vec3(p.x, p.y, p.z);
    float t = dot(v, p3);

    r.w = g * (p.w - t);
    r.x = -v.x*g*p.w + p.x + (gm1*v.x/v2)*v.x*p.x + (gm1*v.x/v2)*v.y*p.y + (gm1*v.x/v2)*v.z*p.z;
    r.y = -v.y*g*p.w + p.y + (gm1*v.y/v2)*v.x*p.x + (gm1*v.y/v2)*v.y*p.y + (gm1*v.y/v2)*v.z*p.z,
    r.z = -v.z*g*p.w + p.z + (gm1*v.z/v2)*v.x*p.x + (gm1*v.z/v2)*v.y*p.y + (gm1*v.z/v2)*v.z*p.z;
*/

/*	STEP 3
    vec3 p3 = vec3(p.x, p.y, p.z);
    float t = dot(v, p3);

    r.w = g * (p.w - t);
    r.x = -v.x*g*p.w + p.x + gm1*v.x/v2 * t;
    r.y = -v.y*g*p.w + p.y + gm1*v.y/v2 * t,
    r.z = -v.z*g*p.w + p.z + gm1*v.z/v2 * t;
*/

/*	STEP 4
    vec3 p3 = vec3(p.x, p.y, p.z);
    float t = dot(v, p3);

    r.w = g * (p.w - t);
    r.x = v.x * (-g*p.w + gm1*t/v2) + p.x;
    r.y = v.y * (-g*p.w + gm1*t/v2) + p.y,
    r.z = v.z * (-g*p.w + gm1*t/v2) + p.z;
*/

/*	STEP 5 
    vec3 p3 = vec3(p.x, p.y, p.z);
    float t = dot(v, p3);
    float t2 = gm1*t/v2 - g*p.w;

    r.w = g * (p.w - t);
    r.x = v.x*t2 + p.x;
    r.y = v.y*t2 + p.y,
    r.z = v.z*t2 + p.z;
*/

/*	STEP 6  */

    vec3 p3 = vec3(p.x, p.y, p.z);
    float t = dot(v, p3);
    float t2 = gm1*t/v2 - g*p.w;
    r = vec4( v*t2 + p3, g * (p.w - t));

    return r;
}