input.x0 = 7; // min=1.01, max=20, step=0.001

// The equations of Lotka-Volterra are a simple model of
// prey-predator population evolution.
// x represents the amount of preys
// y represents the amount of predators
// Their evolution equation is
// x' = x * (1 - y)
// y' = y * (x - 1)
// When there are few preys (x<1), the number of predators decreases.
// Once there are few predators (y<1), the prey population is able to increase exponentially...
// except this allows the number of predators to increase again, and eat them, leading to fewer preys,
// and the cycle repeats...
//
// This specific model has stable cycles: from any starting position (x0,y0) with positive values,
// the trajectory will return to (x0,y0) after one cycle, thus
// the behavior will be periodic (but the period depends on the starting position).
// In particular there is no "damping" toward the fixed point (x=1, y=1)
// 
// If we start from the point (x0, 1), larger values of x0 mean that the cycle will have stronger variations,
// but also a longer period.

// In terms of sound, I think it's rather cool,
// but the fact that it detunes depending on x0 is inconvenient for using it in synth code.


// Retune the synth !
function calc_period_lk(x0, a0) {
    // How many samples does it take to run one period?
    var x = x0;
    var y = 1 + 0.5 * (x-1) * a0;
    var s = 0;
    while(y >= 1) {
        x += (x * (1-y)) * a0;
        y += (y * (x-1)) * a0;
        s++;
    }
    var dy = 0;
    while(y < 1) {
        x += (x * (1-y)) * a0;
        dy = (y * (x-1)) * a0;
        y += dy;
        s++;
    }
    // Sub-sample correction
    var overshoot = (y-1) / dy;
    return s - overshoot; 
}

const lk = synth.def(class {
    constructor(options) {
        this.x0 = input.x0;
        this.calc_a0(options);
        this.x = this.x0;
        this.y = 1 + 0.5 * (this.x-1) * this.a0;
    }
    calc_a0(options) {
        // Change the coefficient a0 so that the frequency is correct despite the trajectory being different
        var freq = midi_to_hz(options.note);
        var a0 = 2*Math.PI * freq * ditty.dt;
        var spl = calc_period_lk(this.x0, a0); // current period duration (in samples)
        debug.log("spl", spl);
        var spl_target = 1 / (freq * ditty.dt); // desired period duration
        debug.log("spl_target", spl_target);
        var factor = (spl / spl_target);
        this.a0 = a0 * factor; // Just multiply the speed by this factor and hope for the best
        debug.log("factor", factor);
        var spl2 = calc_period_lk(this.x0, this.a0); // Now the period is a bit closer, but still not perfect
        debug.log("spl2", spl2);
        var factor2 = (spl2 / spl_target);
        this.a0 *= factor2; // So we correct the factor again
        var spl3 = calc_period_lk(this.x0, this.a0); // Almost perfect now
        debug.log("spl3", spl3);
        var factor3 = (spl3 / spl_target);
        this.a0 *= factor3; // One last time
        
        // It can be nice to use adaptive oversampling for large a0,
        // however this implementation is not compatible with the retuning.
        this.os = 1;
        /*
        while(this.a0 > 0.01) {
            this.os *= 2;
            this.a0 /= 2;
        }
        debug.log("os",this.os);
        */
    }
    process(note,env,tick,options) {
        if(this.x0 != input.x0) { // Reset on input changd
            this.x0 = input.x0;
            this.x = this.x0;
            this.calc_a0(options);
            this.y = 1 + 0.5 * (this.x-1) * this.a0;
        }
        for(let i=0; i<this.os; i++) {
            // Störmer-Verlet integrator 
            // (symplectic numerical integrator, so that it doesn't explode nor dampen)
            this.x += (this.x * (1-this.y)) * this.a0;
            this.y += (this.y * (this.x-1)) * this.a0;
        }
        return this.y - 1;
    }
}, {env:one, amp:0.03});

lk.play(c3);