sclass ContinuousOscillators {

  macro dprint { if (debug) printVars } bool debug;

  int channels = 1;

  double sampleRate = 48000;

  double currentSample;

  IQuerySound sound;

  IAudioSample audio;

  SlidingWindow window;

  bool renderUsingSine; // or square wave

  bool realOnly, imagOnly; // for testing

  double maxRelFreq = 0.25; // maximum oscillator frequency relative to sample rate

  long haarFeaturesChecked;

  persistable class Oscillator extends HasKey<Frequency> {

    double currentPeriodStart;

    // all recorded intensities

    new TreeMap<DoubleRange, Channels<Complex>> intensities;

    double interval; // length of period in samples

    *(Frequency f) {

      super(f);

      interval = frequency().interval()*sampleRate;

    }

    Frequency frequency() { ret key(); }

    double interval() { ret interval; }

    void record_impl(DoubleRange r, Channels<Complex> c) {

      intensities.put(r, c);

    }

    void measure() {

      DoubleRange p1 = doubleRangeWithLength(currentPeriodStart, interval);

      AudioHaarFeature haar = new(audio, p1);

      record(this, p1, haar.getComplex());

      haarFeaturesChecked += 2;

      schuleNextOscillation();

    }

    void initial() {

      completionsForSample.put(currentPeriodStart, this);

    }

    void schuleNextOscillation() {

      completionsForSample.put(currentPeriodStart = phase360(), this);

    }

    double phase360() { ret currentPeriodStart + interval(); }

  } // end of Oscillator

  TreeHasKeyMap<Frequency, Oscillator> oscillators;

  new TreeMultiMap<Double, Oscillator> completionsForSample;

  *() {}

  *(IAudioSample audio, double *currentSample) { setAudio(audio); }

  *(IAudioSample audio) { setAudio(audio); } // start at time 0

  void setFrequencies(Iterable<Frequency> frequencies) {

    oscillators = new TreeHasKeyMap<Frequency, Oscillator>(map(f -> new Oscillator(f), filter(frequencies,

      f -> between(f!, 1, sampleRate*maxRelFreq))));

  }

  void setAudio(IAudioSample audio) {

    this.audio = audio;

    window = optCast SlidingWindow(audio);

  }

  int nFrequencies() { ret oscillators.size(); }

  void startOscillators {

    for (o : oscillators) o.initial();

  }

  void stepTo(DoubleRange timeRange) {

    while ping (true) {

      Double t = firstKey(completionsForSample);

      if (t != null && t <= timeRange.end) {

        var list = completionsForSample.getAndClear(t);

        dprint(+t, +list);

        for (Oscillator o : list)

          o.measure();

      } else

        break;

    }

  }

  void record(Oscillator o, DoubleRange r, Channels<Complex> c) {

    if (o == null) ret;

    o.record_impl(r, c);

  }

  Frequency lowestFrequency() { ret oscillators.firstKey(); }

  Frequency highestFrequency() { ret oscillators.lastKey(); }

  int minWindowSize() {

    ret iceil(sampleRate*lowestFrequency().interval());

  }

  void makeSmallestWindow(IQuerySound sound) {

    this.sound = sound;

    setAudio(new SlidingWindow(channels, sampleRate, sound, 0, minWindowSize());

  }

  void makeFullWindow(IQuerySound sound, IntRange r) {

    this.sound = sound;

    setAudio(new SlidingWindow(channels, sampleRate, sound,

      r.start, r.length()));

  }

  double windowSize() { ret window.length(); }

  S stats() { ret renderVars(+audio, +haarFeaturesChecked); }

  // how much horizontal area have we covered

  L<Double> coverageByFreq() {

    ret mapReversed(oscillators, o -> totalLengthOfDoubleRanges(keys(o.intensities)));

  }

  DoubleRange windowBounds() { ret window.bounds(); }

  IIntegralImage imageColumn(DoubleRange timeRange) {

    int n = nFrequencies();

    double[] col = rawImageColumn(timeRange);

    double[] col2 = normalizeDoubles(col, 255);

    ret bwIntegralImageFromFunction(1, n, (x, y) -> ifloor(col2[y]));

  }

  double[] rawImageColumn(DoubleRange timeRange) {

    int n = nFrequencies();

    var l = reversedList(oscillators);

    double[] col = new[n];

    for y to n: {

      TreeMap<DoubleRange, Channels<Complex>> l2 = l.get(y).intensities;

      // find right period

      DoubleRange period = l2.floorKey(timeRange);

      if (period == null) continue;

      double amplitude = l2.get(period).first().abs();

      col[y] = amplitude;

    }

    ret col;

  }

  MakesBufferedImage imageForWindow(double pixelsPerSecond) {

    new L<IIntegralImage> images;

    double len = windowBounds().length()/sampleRate;

    int pixels = iround(len*pixelsPerSecond);

    for i to pixels: {

      DoubleRange range = doubleRange(

        i/pixelsPerSecond*sampleRate+windowBounds().start(),

        (i+1)/pixelsPerSecond*sampleRate+windowBounds().start());

      //print(+range);

      images.add(imageColumn(range));

    }

    ret mergeBWImagesHorizontally(map toBWImage(images), spacing := 0);

  }

  double[] toAudio(IntRange sampleRange, Iterable<Oscillator> oscillatorsToRender default oscillators) {

    double[] samples = new[l(sampleRange)];

    class PerFreq {

      Oscillator o;

      DoubleRange period;

      Channels<Complex> intensity;

      *(Oscillator *o) {

        setPeriod(o.intensities.floorKey(toDoubleRange(sampleRange)));

        //printVars("PerFreq", f := o.frequency(), +sampleRange, +period, firstPeriod := firstKey(o.intensities));

      }

      void setPeriod(DoubleRange period) {

        this.period = period;

        intensity = period == null ?: o.intensities.get(period);

      }

      void nextPeriod(double t) {

        while (period != null && t >= period.end)

          setPeriod(o.intensities.higherKey(period));

      }

    }

    L<PerFreq> perFreqs = map(oscillatorsToRender, o -> new PerFreq(o));

    for i over samples: {

      double t = sampleRange.start+i; // t is time in samples

      double sum = 0;

      for (pf : perFreqs) {

        pf.nextPeriod(t);

        bool shouldPrint = toAudio_shouldPrint(i, l(samples));

        if (pf.intensity != null) {

          Complex c = div(pf.intensity.first(), l(pf.period));

          double t2 = t-pf.period.start;

          double frac = t2/pf.o.interval;

          //double phiRE = frac*pi()*2;

          //double phiIM = phiRE+pi()/2;

          if (shouldPrint)

            printVars(f := pf.o.frequency(), +i, interval := formatDouble1(pf.o.interval)/*, interval2 := l(pf.period)*/, +t, +t2, frac := formatDouble2(frac),

              c := renderComplexWithAngle(div(c, 32768)));

          if (realOnly) c = complex(c.re());

          else if (imagOnly) c = complex(0, c.im());

          sum += renderFrequencySample(c, frac);

        } else

          if (shouldPrint)

            printVars(f := pf.o.frequency(), +i, interval := pf.o.interval, +t, c := 0);

      }

      samples[i] = sum;

    }

    ret samples;

  }

  double renderFrequencySample(Complex c, double frac) {

    ret (renderUsingSine

      ? new RenderFrequencySample_Sine

      : new RenderFrequencySample_SquareWave).get(c, frac);

  }

  swappable bool toAudio_shouldPrint(int i, int n) { false; }

}