if(GPU.GPU == undefined) gpu = new GPU()
else gpu = new GPU.GPU()
var width = 800; var height = 800;

var paramNames = ['theta', 'dist', 'scale', 'SR', 'discR', 'discW', 'iters']
var paramInputs = paramNames.map((n) => document.getElementById(n))
var autoUpdate = document.getElementById('autoupdate')

var kernel = gpu.createKernel(function(pTheta, pDist, pScale, RS, pDiscR, pDiscW, iters)
{
    // Tunable parameters
    let x = (this.thread.x - 400) / 1 * 0.0001
    let y = (this.thread.y - 400) / 1 * 0.0001
    
    let startPoint = [pDist * Math.cos(pTheta / 180 * 3.141592654), 0, pDist * Math.sin(pTheta / 180 * 3.141592654)];
    let startVelocity = [-1.0, x * pScale + 0.0007, y * pScale];
    let c = Math.cos(-pTheta / 180 * 3.141592654);
    let s = Math.sin(-pTheta / 180 * 3.141592654);
    startVelocity = [startVelocity[0] * c + startVelocity[2] * s, startVelocity[1], startVelocity[2] * c - startVelocity[0] * s]

    // Start of general computation
    let planeNormal = v3normalize(v3cross(startPoint, startVelocity));

    let startDist = v3abs(startPoint);

    let dr = v3dot(startPoint, startVelocity) / startDist;
    let dph = Math.sqrt(v3dot(startVelocity, startVelocity) - dr * dr) / startDist;
    let dt = Math.sqrt(
        (dr * dr / (1.0 - RS / startDist) + startDist * startDist * dph * dph) / (1.0 - RS / startDist)
    )
    let point = [startDist, 0.0];
    let velocity = [dr / dt, dph / dt];
    let finite = true;
    let timestep = 10.0 / iters;

    for(let i = 0; i < iters; i++)
    {
        let acc = getAcceleration(point, velocity, RS);
        point = v2add(point, v2mul(velocity, timestep));
        velocity = v2add(velocity, v2mul(acc, timestep));

        let dratio = (point[0] / startDist)
        let cartesianPoint = v3add(v3mul(startPoint, Math.cos(point[1]) * dratio), v3mul(v3cross(planeNormal, startPoint), Math.sin(point[1]) * dratio))
        if(Math.sqrt(cartesianPoint[0]*cartesianPoint[0] + cartesianPoint[1]*cartesianPoint[1]) < pDiscR && Math.abs(cartesianPoint[2]) < pDiscW)
        {
            let c1 = [236,52,0]
            let c2 = [255,225,67]
            let t = 1.0 - (i / iters)
            let r = v3mul(v3add(c1, v3mul(v3sub(c2, c1), t * 5 - 0.8)), 1 / 255)
            this.color(r[0], r[1], r[2]);
            // this.color(c, c, c)
            finite = false
            break
        }
        if(point[0] < RS)
        {
            this.color(1.0, 0.0, 0.0);
            finite = false;
            break
        }
    }
    if(finite)
        this.color(0.0, 0.0, 0.0);
}, {output: [width, height], debug: true, graphical: true, loopMaxIterations: 1000000, tactice: "highp"}).setFunctions([
    v2add, v2sub, v2mul, 
    v3add, v3sub, v3mul, v3div, v3dot, v3cross, v3abs, v3normalize, 
getAcceleration])

function tryUpdate()
{
    if(autoUpdate.checked) run();
}

function run()
{
    let params = paramInputs.map((el) => parseFloat(el.value))
    kernel(...params);

    try {
    } catch (error) {
        console.log(error);
    }
    
    const canvas = kernel.canvas;
    document.getElementById('here').appendChild(canvas);
}

function getAcceleration(point, velocity, RS)
{
    let r = point[0]

    let Gtrt = RS / (2 * r * (r - RS))

    let Grrr = -Gtrt
    let Grtt = RS * (r - RS) / (2 * r*r*r)
    let Grphph = (RS - r)

    let Gphrph = 1 / r

    let vr = velocity[0]
    let vph = velocity[1]

    return [
        - Grrr * vr * vr - Grtt - Grphph * vph * vph + 2 * Gtrt * vr * vr,
        - 2 * Gphrph * vr * vph + 2 * Gtrt * vr * vph
    ]
}

function v2add(a, b)
{
    return [a[0] + b[0], a[1] + b[1]];
}

function v2sub(a, b)
{
    return [a[0] - b[0], a[1] - b[1]];
}

function v2mul(a, s)
{
    return [a[0] * s, a[1] * s];
}

function v3add(a, b)
{
    return [a[0] + b[0], a[1] + b[1], a[2] + b[2]];
}

function v3sub(a, b)
{
    return [a[0] - b[0], a[1] - b[1], a[2] - b[2]];
}

function v3mul(a, s)
{
    return [a[0] * s, a[1] * s, a[2] * s];
}

function v3div(a, s)
{
    return [a[0] / s, a[1] / s, a[2] / s];
}

function v3dot(a, b)
{
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

function v3cross(a, b)
{
    return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
}

function v3abs(x)
{
    // return Math.hypot(x[0], x[1], x[2]);
    return Math.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
}

function v3normalize(vec)
{
    return v3div(vec, v3abs(vec));
}