#ifndef __OSCILLATOR_H__
#define __OSCILLATOR_H__

/* PolyBLEP oscillator inspired by: https://www.metafunction.co.uk/post/all-about-digital-oscillators-part-2-blits-bleps
 * Further info about BLEP at: https://pbat.ch/sndkit/blep/ 
 */

#include <cstdlib>
#include <math.h>

#include "util.h"

class Oscillator {
public:
    enum Mode { MODE_SINE = 0, MODE_SAW, MODE_SQUARE };

    Mode const * mode;
    float phase;            // current waveform phase angle (radians)
    float value;            // current amplitude value
    float driftAmount;
    float driftValue;

    Oscillator() : phase(randPhase()), value(0), driftAmount(0.001f), driftValue(0) {}

    void assign(Mode const * mode);

    // Generate next output sample and advance the phase angle
    inline float tick(float cv) {
        const float phaseStep = (float) (2 * CV_FREQ_MIN) * exp2f(cv) * (float) SAMPLE_RATE_INV;
        
        const float t = phase * (float) PIx2_INV; // Define half phase
        const float dt = phaseStep * (float) PIx2_INV;
        
        if (*mode == MODE_SINE) {
            value = sinf(phase); // No harmonics in sine so no aliasing!! No Poly BLEPs needed!
        } else if (*mode == MODE_SAW) {
            value = (2.0f * t) - 1.0f; // Render naive waveshape
            value -= polyBlep(t, dt); // Layer output of Poly BLEP on top
        } else if (*mode == MODE_SQUARE) {
            if (phase < (float) M_PI) {     
                value = 1.0; // Flip 
            } else { 
                value = -1.0; // Flop
            }
            value += polyBlep(t, dt); // Layer output of Poly BLEP on top (flip)
            value -= polyBlep(fmodf(t + 0.5f, 1.0f), dt); // Layer output of Poly BLEP on top (flop)
        }

        driftValue += 0.005f * whiteNoise() - 0.00001f * driftValue;
        phase += phaseStep * (1.0f + driftAmount * driftValue);

        if(phase >= (float) PIx2) { // wrap if phase angle >=360º
            phase -= PIx2;
        }
        
        return value;
    }

    inline float polyBlep(float t, float dt) {
        // t-t^2/2 +1/2
        // 0 < t <= 1
        // discontinuities between 0 & 1
        if (t < dt) {
            t /= dt;
            return t + t - t * t - 1.0f;
        }

        // t^2/2 +t +1/2
        // -1 <= t <= 0
        // discontinuities between -1 & 0
        else if (t > 1.0f - dt) {
            t = (t - 1.0f) / dt;
            return t * t + t + t + 1.0f;
        }

        // no discontinuities 
        // 0 otherwise
        else return 0.0;
    }
};

#endif