//############################## oscillators.lib ######################################
// This library contains a collection of sound generators. Its official prefix is `os`.
//########################################################################################

/************************************************************************
************************************************************************
FAUST library file, GRAME section

Except where noted otherwise, Copyright (C) 2003-2017 by GRAME,
Centre National de Creation Musicale.
----------------------------------------------------------------------
GRAME LICENSE

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2.1 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA.

EXCEPTION TO THE LGPL LICENSE : As a special exception, you may create a
larger FAUST program which directly or indirectly imports this library
file and still distribute the compiled code generated by the FAUST
compiler, or a modified version of this compiled code, under your own
copyright and license. This EXCEPTION TO THE LGPL LICENSE explicitly
grants you the right to freely choose the license for the resulting
compiled code. In particular the resulting compiled code has no obligation
to be LGPL or GPL. For example you are free to choose a commercial or
closed source license or any other license if you decide so.
************************************************************************
************************************************************************/

declare name "Faust Oscillator Library";
declare version "0.1";

ma = library("maths.lib");
ba = library("basics.lib");
fi = library("filters.lib");

// This library contains platform specific constants 
pl = library("platform.lib");

//=========================Wave-Table-Based Oscillators===================================
//========================================================================================


//-----------------------`(os.)sinwaveform`------------------------
// Sine waveform ready to use with a `rdtable`.
//
// #### Usage
//
// ```
// sinwaveform(tablesize) : _
// ```
//
// Where:
//
// * `tablesize`: the table size
//------------------------------------------------------------
sinwaveform(tablesize) = float(ba.time)*(2.0*ma.PI)/float(tablesize) : sin;

//-----------------------`(os.)coswaveform`------------------------
// Cosine waveform ready to use with a `rdtable`.
//
// #### Usage
//
// ```
// coswaveform(tablesize) : _
// ```
//
// Where:
//
// * `tablesize`: the table size
//------------------------------------------------------------
coswaveform(tablesize) = float(ba.time)*(2.0*ma.PI)/float(tablesize) : cos;

//-----------------------`(os.)phasor`------------------------
// A simple phasor to be used with a `rdtable`.
// `phasor` is a standard Faust function.
//
// #### Usage
//
// ```
// phasor(tablesize,freq) : _
// ```
//
// Where:
//
// * `tablesize`: the table size
// * `freq`: the frequency of the phasor (Hz)
//------------------------------------------------------------
phasor(tablesize,freq) = freq/float(ma.SR) : (+ : ma.decimal) ~ _ : *(float(tablesize));

//-----------------------`(os.)hs_phasor`------------------------
// Hardsyncing phasor to be used with an `rdtable`.
//
// #### Usage
//
// ```
// hs_phasor(tablesize,freq,c) :  _
// ```
//
// Where:
//
// * `tablesize`: the table size
// * `freq`: the frequency of the phasor (Hz)
// * `c`: a clock signal, `c>0` resets phase to 0
//---------------------------------------------------------
// Author: Mike Olsen
hs_phasor(tablesize,freq,c) = inc : (+ : d)~(-(_<:(_,*(_,clk)))) : *(tablesize)
with {
	clk = c>0;
	d   = ma.decimal;
	inc = freq/float(ma.SR);
};

//-----------------------`(os.)hsp_phasor`------------------------
// Hardsyncing phasor with selectable phase to be used with an `rdtable`.
//
// #### Usage
//
// ```
// hsp_phasor(tablesize,freq,p,c)
// ```
//
// Where:
//
// * `tablesize`: the table size
// * `freq`: the frequency of the phasor (Hz)
// * `p`: phase init (0 <= `p` <=1)
// * `c`: a clock signal, `c>0` resets phase to 0
//---------------------------------------------------------
// Author: Christophe Lebreton
hsp_phasor(tablesize,freq,p,c) = inc : (+ : d)~(-(_<:(_,*(_,clk))):+(clk*p)) : *(tablesize)
with {
    clk = c>0;
    d   = ma.decimal;
    inc = freq/float(ma.SR);
};

//-----------------------`(os.)oscsin`------------------------
// Sine wave oscillator.
// `oscsin` is a standard Faust function.
//
// #### Usage
//
// ```
// oscsin(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency of the wave (Hz)
//------------------------------------------------------------
oscsin(freq) = rdtable(tablesize, sinwaveform(tablesize), int(phasor(tablesize,freq)))
with {
	tablesize = pl.tablesize;
};

//-----------------------`(os.)hs_oscsin`------------------------
// Sin lookup table with hardsyncing phase.
//
// #### Usage
//
// ```
// hs_oscsin(freq,c) : _
// ```
//
// Where:
//
// * `freq`: the fundamental frequency of the phasor
// * `c`: a clock signal, `c>0` resets phase to 0
//---------------------------------------------------------
// Author: Mike Olsen
hs_oscsin(freq,c) = rdtable(tablesize, sinwaveform(tablesize), int(hs_phasor(tablesize,freq,c)))
with {
	tablesize = pl.tablesize;
};

//-----------------------`(os.)osccos`------------------------
// Cosine wave oscillator.
//
// #### Usage
//
// ```
// osccos(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency of the wave (Hz)
//------------------------------------------------------------
osccos(freq) = rdtable(tablesize, coswaveform(tablesize), int(phasor(tablesize,freq)))
with {
	tablesize = pl.tablesize;
};

//-----------------------`(os.)oscp`------------------------
// A sine wave generator with controllable phase.
//
// #### Usage
//
// ```
// oscp(freq,phase) : _
// ```
//
// Where:
//
// * `freq`: the frequency of the wave (Hz)
// * `phase`: the phase in radian
//------------------------------------------------------------
oscp(freq,phase) = oscsin(freq) * cos(phase) + osccos(freq) * sin(phase);

//-----------------------`(os.)osci`------------------------
// Interpolated phase sine wave oscillator.
//
// #### Usage
//
// ```
// osci(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency of the wave (Hz)
//------------------------------------------------------------
osci(freq) = s1 + d * (s2 - s1)
with {
	tablesize = pl.tablesize;
	i = int(phasor(tablesize,freq));
	d = ma.decimal(phasor(tablesize,freq));
	s1 = rdtable(tablesize+1,sinwaveform(tablesize),i);
	s2 = rdtable(tablesize+1,sinwaveform(tablesize),i+1);
};

// end GRAME section
//########################################################################################
/************************************************************************
FAUST library file, jos section

Except where noted otherwise, The Faust functions below in this
section are Copyright (C) 2003-2017 by Julius O. Smith III <jos@ccrma.stanford.edu>
([jos](http://ccrma.stanford.edu/~jos/)), and released under the
(MIT-style) [STK-4.3](#stk-4.3-license) license.

The MarkDown comments in this section are Copyright 2016-2017 by Romain
Michon and Julius O. Smith III, and are released under the
[CCA4I](https://creativecommons.org/licenses/by/4.0/) license (TODO: if/when Romain agrees)

************************************************************************/

//===============================LFOs===============================
// Low-Frequency Oscillators (LFOs) have prefix `lf_`
// (no aliasing suppression, which is not audible at LF).
//==================================================================

//--------`(os.)lf_imptrain`----------
// Unit-amplitude low-frequency impulse train.
// `lf_imptrain` is a standard Faust function.
//
// #### Usage
//
// ```
// lf_imptrain(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_imptrain(freq) = lf_sawpos(freq)<:-(mem)<0; // definition below

//--------`(os.)lf_pulsetrainpos`----------
// Unit-amplitude nonnegative LF pulse train, duty cycle between 0 and 1.
//
//
// #### Usage
//
// ```
// lf_pulsetrainpos(freq,duty) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `duty`: duty cycle between 0 and 1
//------------------------------------------------------------
lf_pulsetrainpos(freq,duty) = float(lf_sawpos(freq) <= duty);

//pulsetrainpos = lf_pulsetrainpos; // for backward compatibility

//--------`(os.)lf_pulsetrain`----------
// Unit-amplitude zero-mean LF pulse train, duty cycle between 0 and 1.
//
// #### Usage
//
// ```
// lf_pulsetrain(freq,duty) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
// * `duty`: duty cycle between 0 and 1
//------------------------------------------------------------
lf_pulsetrain(freq,duty) = 2.0*lf_pulsetrainpos(freq,duty) - 1.0;

//--------`(os.)lf_squarewavepos`----------
// Positive LF square wave in [0,1]
//
// #### Usage
//
// ```
// lf_squarewavepos(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_squarewavepos(freq) = lf_pulsetrainpos(freq,0.5);
// squarewavepos = lf_squarewavepos; // for backward compatibility


//--------`(os.)lf_squarewave`----------
// Zero-mean unit-amplitude LF square wave.
// `lf_squarewave` is a standard Faust function.
//
// #### Usage
//
// ```
// lf_squarewave(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_squarewave(freq) = 2*lf_squarewavepos(freq) - 1;
// squarewave = lf_squarewave; // for backward compatibility


//--------`(os.)lf_trianglepos`----------
// Positive unit-amplitude LF positive triangle wave.
//
// #### Usage
//
// ```
// lf_trianglepos(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
lf_trianglepos(freq) = 1-abs(saw1(freq)); // saw1 defined below

//----------`(os.)lf_triangle`----------
// Positive unit-amplitude LF triangle wave
// `lf_triangle` is a standard Faust function.
//
// #### Usage
//
// ```
// lf_triangle(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency in Hz
//------------------------------------------------------------
// Author: Bart Brouns
// License: STK-4.3
lf_triangle(freq) = 2*lf_trianglepos(freq) - 1;

//================== Low Frequency Sawtooths ====================
// Sawtooth waveform oscillators for virtual analog synthesis et al.
// The 'simple' versions (`lf_rawsaw`, `lf_sawpos` and `saw1`), are mere samplings of
// the ideal continuous-time ("analog") waveforms.  While simple, the
// aliasing due to sampling is quite audible.  The differentiated
// polynomial waveform family (`saw2`, `sawN`, and derived functions)
// do some extra processing to suppress aliasing (not audible for
// very low fundamental frequencies).  According to Lehtonen et al.
// (JASA 2012), the aliasing of `saw2` should be inaudible at fundamental
// frequencies below 2 kHz or so, for a 44.1 kHz sampling rate and 60 dB SPL
// presentation level;  fundamentals 415 and below required no aliasing
// suppression (i.e., `saw1` is ok).
//=====================================================================

//-----------------`(os.)lf_rawsaw`--------------------
// Simple sawtooth waveform oscillator between 0 and period in samples.
//
// #### Usage
//
// ```
// lf_rawsaw(periodsamps)
// ```
//
// Where:
//
// * `periodsamps`: number of periods per samples
//---------------------------------------------------------
lf_rawsaw(periodsamps) = (_,periodsamps : fmod) ~ +(1.0);

//-----------------`(os.)lf_sawpos_phase`--------------------
// Simple sawtooth waveform oscillator between 0 and 1
// with phase control.
//
// #### Usage
//
// ```
// lf_sawpos_phase(phase,freq)
// ```
//
// Where:
//
// * `phase`: phase
// * `freq`: frequency
//---------------------------------------------------------
lf_sawpos_phase(phase,freq) = (+(phase-phase') : ma.frac) ~ +(freq/ma.SR);

//-----------------`(os.)lf_sawpos`--------------------
// Simple sawtooth waveform oscillator between 0 and 1.
//
// #### Usage
//
// ```
// lf_sawpos(freq)
// ```
//
// Where:
//
// * `freq`: frequency
//
//---------------------------------------------------------
// Author: Bart Brouns
// License: STK-4.3
// MarkDown: Romain Michon
lf_sawpos(freq) = ma.frac ~ +(freq/ma.SR);

//-----------------`(os.)lf_sawpos_reset`--------------------
// Simple sawtooth waveform oscillator between 0 and 1.
// with reset.
//
// #### Usage
//
// ```
// lf_sawpos_reset(freq,reset)
// ```
//
// Where:
//
// * `freq`: frequency
// * `reset`: reset the oscillator to 0
//
//---------------------------------------------------------
// Author: Bart Brouns
// License: STK-4.3
lf_sawpos_reset(freq,reset) = ma.frac * (reset == 0) ~ +(freq/ma.SR);

//-----------------`(os.)lf_sawpos_phase_reset`--------------------
// Simple sawtooth waveform oscillator between 0 and 1.
// with phase control and reset.
//
// #### Usage
//
// ```
// lf_sawpos_phase_reset(freq,phase,reset)
// ```
//
// Where:
//
// * `freq`: frequency
// * `phase`: phase between 0 and 1
// * `reset`: reset the oscillator to 0
//
//---------------------------------------------------------
// Author: Bart Brouns
// License: STK-4.3
lf_sawpos_phase_reset(freq,phase,reset) = lf_sawpos_reset(freq,reset) + phase : ma.frac;

//-----------------`(os.)lf_saw`--------------------
// Simple sawtooth waveform.
// `lf_saw` is a standard Faust function.
//
// #### Usage
//
// ```
// lf_saw(freq)
// ```
//
// Where:
//
// * `freq`: frequency
//---------------------------------------------------------
// Author: Bart Brouns
// License: STK-4.3
saw1(freq) = 2.0 * lf_sawpos(freq) - 1.0;
lf_saw(freq) = saw1(freq);

//================== Bandlimited Sawtooth ====================
//-----------------`(os.)sawN`--------------------
// Bandlimited Sawtooth
//
// `sawN(N,freq)`, `sawNp`, `saw2dpw(freq)`, `saw2(freq)`, `saw3(freq)`,
// `saw4(freq)`, `saw5(freq)`, `saw6(freq)`, `sawtooth(freq)`, `saw2f2(freq)`
// `saw2f4(freq)`
//
// #### Method 1 (`saw2`)
//
// Polynomial Transition Regions (PTR) (for aliasing suppression).
//
// ##### References
//
// * Kleimola, J.; Valimaki, V., "Reducing Aliasing from Synthetic Audio
// 	Signals Using Polynomial Transition Regions," in Signal Processing
// 	Letters, IEEE , vol.19, no.2, pp.67-70, Feb. 2012
// * <https://aaltodoc.aalto.fi/bitstream/handle/123456789/7747/publication6.pdf?sequence=9>
// * <http://research.spa.aalto.fi/publications/papers/spl-ptr/>
//
// #### Method 2 (`sawN`)
//
// Differentiated Polynomial Waves (DPW) (for aliasing suppression).
//
// ##### Reference
//
// "Alias-Suppressed Oscillators based on Differentiated Polynomial Waveforms",
// Vesa Valimaki, Juhan Nam, Julius Smith, and Jonathan Abel,
// IEEE Tr. Acoustics, Speech, and Language Processing (IEEE-ASLP),
// Vol. 18, no. 5, May 2010.
//
// #### Other Cases
//
// Correction-filtered versions of `saw2`: `saw2f2`, `saw2f4`
// The correction filter compensates "droop" near half the sampling rate.
// See reference for sawN.
//
// #### Usage
//
// ```
// sawN(N,freq) : _
// sawNp(N,freq,phase) : _
// saw2dpw(freq) : _
// saw2(freq) : _
// saw3(freq) : _ // based on sawN
// saw4(freq) : _ // based on sawN
// saw5(freq) : _ // based on sawN
// saw6(freq) : _ // based on sawN
// sawtooth(freq) : _ // = saw2
// saw2f2(freq) : _
// saw2f4(freq) : _
// ```
//
// Where:
//
// * `N`: polynomial order
// * `freq`: frequency in Hz
// * `phase`: phase
//===================================================================
// --- sawN for N = 1 to 6 ---
//We can do 6, but 5 and 6 have noise at low fundamentals: MAX_SAW_ORDER = 6; MAX_SAW_ORDER_NEXTPOW2 = 8;
MAX_SAW_ORDER = 4; MAX_SAW_ORDER_NEXTPOW2 = 8; // par cannot handle the case of 0 elements
sawN(N,freq) = saw1l : poly(Nc) : D(Nc-1) : gate(Nc-1)
with {
  Nc = max(1,min(N,MAX_SAW_ORDER));
  clippedFreq = max(20.0,abs(freq)); // use lf_sawpos(freq) for LFOs (freq < 20 Hz)
  saw1l = 2*lf_sawpos(clippedFreq) - 1; // zero-mean, amplitude +/- 1
  // Also note the availability of lf_sawpos_phase above.
  poly(1,x) = x;
  poly(2,x) = x*x;
  poly(3,x) = x*x*x - x;
  poly(4,x) = x*x*(x*x - 2.0);
  poly(5,x) = x*(7.0/3 + x*x*(-10.0/3.0 + x*x));
  poly(6,x) = x*x*(7.0 + x*x*(-5.0 + x*x));
  p0n = float(ma.SR)/clippedFreq; // period in samples
  diff1(x) = (x - x')/(2.0/p0n);
  diff(N) = seq(n,N,diff1); // N diff1s in series
  factorial(0) = 1;
  factorial(i) = i * factorial(i-1);
  D(0) = _;
  D(i) = diff(i)/factorial(i+1);
  gate(N) = *(1@(N)); // delayed step for blanking startup glitch
};

//------------------`(os.)sawNp`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
// --- sawNp for N = 1 to 6 ---
// Phase offset = delay (max 8191 samples is more than one period of audio):
sawNp(N,freq,phase) = sawN(N,freq) : @(max(0,min(8191,int(phase*ma.SR/freq))));

// Special named cases:

//------------------`(os.)saw2dpw`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
// --- sawN ---
saw2dpw(freq) = saw1(freq) <: * <: -(mem) : *(0.25'*ma.SR/freq); // inferior to saw2 below

//------------------`(os.)saw3`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
saw3 = sawN(3); saw4 = sawN(4); saw5 = sawN(5); saw6 = sawN(6);

//------------------`(os.)sawtooth`--------------------------------
// Alias-free sawtooth wave. 2nd order interpolation (based
// on `saw2`).
// `sawtooth` is a standard Faust function.
//
// #### Usage
//
// ```
// sawtooth(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//--------------------------------------------------------
saw2(freq) = y with { // newer PTR version (stateless - freq can vary at any speed)
  p0 = float(ma.SR)/float(max(ma.EPSILON,abs(freq))); // period in samples
  t0 = 1.0/p0; // phase increment
  p = ((_<:(-(1)<:_,_),_) <: selector1,selector2) ~(+(t0)):!,_;
  selector1 = select2(<(0)); // for feedback
  selector2 = select2(<(0), (_<:_,(*(1-p0):+(1)):+), _); // for output
  y = 2*p-1;
};
// --- sawtooth ---
sawtooth = saw2; // default choice for sawtooth signal - see also sawN

//------------------`(os.)saw2f2`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
// --- Correction-filtered versions of saw2: saw2f2, saw2f4 ----
// The correction filter compensates "droop" near half the sampling rate.
// See reference for sawN.
saw2f2 = saw2 : cf2 with {
  cf2 = fi.tf2(1.155704605878911, 0.745184288225518,0.040305967265900,
        0.823765146386639, 0.117420665547108);
};

//------------------`(os.)saw2f4`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
saw2f4 = saw2 : cf4 with {
  cf4 = fi.iir((1.155727435125014, 2.285861038554662,
        1.430915027294021, 0.290713280893317, 0.008306401748854),
        (2.156834679164532, 1.559532244409321, 0.423036498118354,
        0.032080681130972));
};

//=========Bandlimited Pulse, Square, and Impulse Trains============
// Bandlimited Pulse, Square, and Impulse Trains.
//
// `pulsetrainN`, `pulsetrain`, `squareN`, `square`, `imptrain`, `imptrainN`,
// `triangle`, `triangleN`
//
// All are zero-mean and meant to oscillate in the audio frequency range.
// Use simpler sample-rounded lf_* versions above for LFOs.
//
// #### Usage
//
// ```
// pulsetrainN(N,freq,duty) : _
// pulsetrain(freq, duty) : _ // = pulsetrainN(2)
// squareN(N, freq) : _
// square : _ // = squareN(2)
// imptrainN(N,freq) : _
// imptrain : _ // = imptrainN(2)
// triangleN(N,freq) : _
// triangle : _ // = triangleN(2)
// ```
//
// Where:
//
// * `N`: polynomial order
// * `freq`: frequency in Hz
//====================================================================

//------------------`(os.)pulsetrainN`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
pulsetrainN(N,freq,duty) = diffdel(sawN(N,freqC),del) with {
 // non-interpolated-delay version: diffdel(x,del) = x - x@int(del+0.5);
 // linearly interpolated delay version (sounds good to me):
    diffdel(x,del) = x-x@int(del)*(1-ma.frac(del))-x@(int(del)+1)*ma.frac(del);
 // Third-order Lagrange interpolated-delay version (see filters.lib):
 // diffdel(x,del) = x - fdelay3(DELPWR2,max(1,min(DELPWR2-2,ddel)));
 DELPWR2 = 2048; // Needs to be a power of 2 when fdelay*() used above.
 delmax = DELPWR2-1; // arbitrary upper limit on diff delay (duty=0.5)
 SRmax = 96000.0; // assumed upper limit on sampling rate
 fmin = SRmax / float(2.0*delmax); // 23.4 Hz (audio freqs only)
 freqC = max(freq,fmin); // clip frequency at lower limit
 period = (float(ma.SR) / freqC); // actual period
 ddel = duty * period; // desired delay
 del = max(0,min(delmax,ddel));
};

//------------------`(os.)pulsetrain`--------------------------------
// Bandlimited pulse train oscillator. Based on `pulsetrainN(2)`.
// `pulsetrain` is a standard Faust function.
//
// #### Usage
//
// ```
// pulsetrain(freq, duty) : _
// ```
//
// Where:
//
// * `freq`: frequency
// * `duty`: duty cycle between 0 and 1
//--------------------------------------------------------
pulsetrain = pulsetrainN(2);

//------------------`(os.)squareN`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
squareN(N,freq) = pulsetrainN(N,freq,0.5);

//------------------`(os.)square`--------------------------------
// Bandlimited square wave oscillator. Based on `squareN(2)`.
// `square` is a standard Faust function.
//
// #### Usage
//
// ```
// square(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//--------------------------------------------------------
square = squareN(2);

//------------------`(os.)impulse`--------------------------------
// One-time impulse generated when the Faust process is started.
// `impulse` is a standard Faust function.
//
// #### Usage
//
// ```
// impulse : _
// ```
//--------------------------------------------------------
impulse = 1-1';

//------------------`(os.)imptrainN`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
imptrainN(N,freq) = impulse + 0.5*ma.diffn(sawN(N,freq));


//------------------`(os.)imptrain`--------------------------------
// Bandlimited impulse train generator. Based on `imptrainN(2)`.
// `imptrain` is a standard Faust function.
//
// #### Usage
//
// ```
// imptrain(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//--------------------------------------------------------
imptrain = imptrainN(2); // default based on saw2

//------------------`(os.)triangleN`--------------------------------
// TODO: MarkDown doc in comments
//--------------------------------------------------------
triangleN(N,freq) = squareN(N,freq) : fi.pole(p) : *(gain) with {
  gain = 4.0*freq/ma.SR; // for aproximate unit peak amplitude
  p = 0.999;
};

//------------------`(os.)triangle`--------------------------------
// Bandlimited triangle wave oscillator. Based on `triangleN(2)`.
// `triangle` is a standard Faust function.
//
// #### Usage
//
// ```
// triangle(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//--------------------------------------------------------
triangle = triangleN(2); // default based on saw2

//===============================Filter-Based Oscillators=================================
// Filter-Based Oscillators
//
// #### Usage
//
// ```
// osc[b|r|rs|rc|s|w](f), where f = frequency in Hz.
// ```
//
// #### References
//
// * <http://lac.linuxaudio.org/2012/download/lac12-slides-jos.pdf>
// * <https://ccrma.stanford.edu/~jos/pdf/lac12-paper-jos.pdf>
//========================================================================================

//--------------------------`(os.)oscb`--------------------------------
// Sinusoidal oscillator based on the biquad.
//
// #### Usage
//
// ```
// oscb(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//------------------------------------------------------------
oscb(f) = impulse : fi.tf2(1,0,0,a1,1)
with {
  a1 = -2*cos(2*ma.PI*f/ma.SR);
};


//--------------------------`(os.)oscrq`---------------------------
// Sinusoidal (sine and cosine) oscillator based on 2D vector rotation,
//  = undamped "coupled-form" resonator
//  = lossless 2nd-order normalized ladder filter.
//
// #### Usage
//
// ```
// oscrq(freq) : _,_
// ```
//
// Where:
//
// * `freq`: frequency
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
oscrq(f) = impulse : fi.nlf2(f,1); // sine and cosine outputs

//--------------------------`(os.)oscrs`---------------------------
// Sinusoidal (sine) oscillator based on 2D vector rotation,
//  = undamped "coupled-form" resonator
//  = lossless 2nd-order normalized ladder filter.
//
// #### Usage
//
// ```
// oscrs(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
oscrs(f) = impulse : fi.nlf2(f,1) : _,!; // sine

//--------------------------`(os.)oscrc`---------------------------
// Sinusoidal (cosine) oscillator based on 2D vector rotation,
//  = undamped "coupled-form" resonator
//  = lossless 2nd-order normalized ladder filter.
//
// #### Usage
//
// ```
// oscrc(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Normalized_Scattering_Junctions.html>
//------------------------------------------------------------
oscrc(f) = impulse : fi.nlf2(f,1) : !,_; // cosine

oscrp(f,p) = oscrq(f) : *(cos(p)), *(sin(p)) : + ; // p=0 for sine, p=PI/2 for cosine, etc.
oscr = oscrs; // default = sine (starts without a pop)

//--------------------------`(os.)oscs`--------------------------------
// Sinusoidal oscillator based on the state variable filter
// = undamped "modified-coupled-form" resonator
// = "magic circle" algorithm used in graphics.
//------------------------------------------------------------
oscs(f) = (*(-1) : sint(wn) : sintp(wn,impulse)) ~ _
with {
  wn = 2*ma.PI*f/ma.SR; // approximate
  // wn = 2*sin(PI*f/SR); // exact
  sint(x) = *(x) : + ~ _ ; // frequency-scaled integrator
  sintp(x,y) = *(x) : +(y): + ~ _ ; // same + state input
};

//-----------------------`(os.)osc`------------------------
// Default sine wave oscillator (same as [oscsin](#oscsin)).
// `osc` is a standard Faust function.
//
// #### Usage
//
// ```
// osc(freq) : _
// ```
//
// Where:
//
// * `freq`: the frequency of the wave (Hz)
//------------------------------------------------------------
osc	= oscsin;

//================ Waveguide-Resonator-Based Oscillators ================
// Sinusoidal oscillator based on the waveguide resonator `wgr`.
//=======================================================================

//-----------------`(os.)oscw`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`. Unit-amplitude
// cosine oscillator.
//
// #### Usage
//
// ```
// oscwc(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscwc(fr) = impulse : fi.wgr(fr,1) : _,!; // cosine (cheapest at 1 mpy/sample)

//-----------------`(os.)oscws`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`. Unit-amplitude
// sine oscillator.
//
// #### Usage
//
// ```
// oscws(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscws(fr) = impulse : fi.wgr(fr,1) : !,_; // sine (needs a 2nd scaling mpy)

//-----------------`(os.)oscwq`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`.
// Unit-amplitude cosine and sine (quadrature) oscillator.
//
// #### Usage
//
// ```
// oscwq(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscq(fr) = impulse : fi.wgr(fr,1);       // phase quadrature outputs

//-----------------`(os.)oscw`--------------------
// Sinusoidal oscillator based on the waveguide resonator `wgr`.
// Unit-amplitude cosine oscillator (default).
//
// #### Usage
//
// ```
// oscw(freq) : _
// ```
//
// Where:
//
// * `freq`: frequency
//
// #### Reference
//
// * <https://ccrma.stanford.edu/~jos/pasp/Digital_Waveguide_Oscillator.html>
//------------------------------------------------------------
oscw = oscwc;

// end jos section
//########################################################################################
/************************************************************************
FAUST library file, further contributions section

All contributions below should indicate both the contributor and terms
of license.  If no such indication is found, "git blame" will say who
last edited each line, and that person can be emailed to inquire about
license disposition, if their license choice is not already indicated
elsewhere among the libraries.  It is expected that all software will be
released under LGPL, STK-4.3, MIT, BSD, or a similar FOSS license.
************************************************************************/

//===================== Casio CZ Oscillators ==========================
// Oscillators that mimic some of the Casio CZ oscillators.
//
// There are two sets:
// - A set with an index parameter
// - A set with a res parameter
//
// The "index oscillators" outputs a sine wave at index=0 and gets brighter with a higher index.
// There are two versions of the "index oscillators":
// - with P appended to the name: is phase aligned with 'fund:sin'
// - without P appended to the name: has the phase of the original CZ oscillators
//
// The "res oscillators" have a resonant frequency.
// "res" is the frequency of resonance as a factor of the fundamental pitch.
//=====================================================================

//----------`(os.)CZsaw`----------
// Oscillator that mimics the Casio CZ saw oscillator
// `CZsaw` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsaw(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = saw-wave
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsaw(fund, index) = CZ.sawChooseP(fund, index, 0);

//----------`(os.)CZsawP`----------
// Oscillator that mimics the Casio CZ saw oscillator,
// with it's phase aligned to `fund:sin`.
// `CZsawP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsawP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = saw-wave
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsawP(fund, index) = CZ.sawChooseP(fund, index, 1);

//----------`(os.)CZsquare`----------
// Oscillator that mimics the Casio CZ square oscillator
// `CZsquare` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsquare(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = square-wave
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsquare(fund, index) = CZ.squareChooseP(fund, index, 0);

//----------`(os.)CZsquareP`----------
// Oscillator that mimics the Casio CZ square oscillator,
// with it's phase aligned to `fund:sin`.
// `CZsquareP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsquareP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 to 1. 0 = sine-wave, 1 = square-wave
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsquareP(fund, index) = CZ.squareChooseP(fund, index, 1);

//----------`(os.)CZpulse`----------
// Oscillator that mimics the Casio CZ pulse oscillator
// `CZpulse` is a standard Faust function.
//
// #### Usage
//
// ```
// CZpulse(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is closer to a pulse
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZpulse(fund, index) = CZ.pulseChooseP(fund, index, 0);

//----------`(os.)CZpulseP`----------
// Oscillator that mimics the Casio CZ pulse oscillator,
// with it's phase aligned to `fund:sin`.
// `CZpulseP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZpulseP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is closer to a pulse
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZpulseP(fund, index) = CZ.pulseChooseP(fund, index, 1);

//----------`(os.)CZsinePulse`----------
// Oscillator that mimics the Casio CZ sine/pulse oscillator
// `CZsinePulse` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsinePulse(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is a sine minus a pulse
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsinePulse(fund, index) = CZ.sinePulseChooseP(fund, index, 0);

//----------`(os.)CZsinePulseP`----------
// Oscillator that mimics the Casio CZ sine/pulse oscillator,
// with it's phase aligned to `fund:sin`.
// `CZsinePulseP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZsinePulseP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is a sine minus a pulse
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZsinePulseP(fund, index) = CZ.sinePulseChooseP(fund, index, 1);

//----------`(os.)CZhalfSine`----------
// Oscillator that mimics the Casio CZ half sine oscillator
// `CZhalfSine` is a standard Faust function.
//
// #### Usage
//
// ```
// CZhalfSine(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is somewhere between a saw and a square
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZhalfSine(fund, index) = CZ.halfSineChooseP(fund, index, 0);

//----------`(os.)CZhalfSineP`----------
// Oscillator that mimics the Casio CZ half sine oscillator,
// with it's phase aligned to `fund:sin`.
// `CZhalfSineP` is a standard Faust function.
//
// #### Usage
//
// ```
// CZhalfSineP(fund,index) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `index`: the brightness of the oscillator, 0 gives a sine-wave, 1 is somewhere between a saw and a square
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZhalfSineP(fund, index) = CZ.halfSineChooseP(fund, index, 1);

//----------`(os.)CZresSaw`----------
// Oscillator that mimics the Casio CZ resonant saw-tooth oscillator
// `CZresSaw` is a standard Faust function.
//
// #### Usage
//
// ```
// CZresSaw(fund,res) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `res`: the frequency of resonance as a factor of the fundamental pitch.
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZresSaw(fund,res) = CZ.resSaw(fund,res);

//----------`(os.)CZresTriangle`----------
// Oscillator that mimics the Casio CZ resonant triangle oscillator
// `CZresTriangle` is a standard Faust function.
//
// #### Usage
//
// ```
// CZresTriangle(fund,res) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `res`: the frequency of resonance as a factor of the fundamental pitch.
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZresTriangle(fund,res) = CZ.resTriangle(fund,res);

//----------`(os.)CZresTrap`----------
// Oscillator that mimics the Casio CZ resonant trapeze oscillator
// `CZresTrap` is a standard Faust function.
//
// #### Usage
//
// ```
// CZresTrap(fund,res) : _
// ```
//
// Where:
//
// * `fund`: a saw-tooth waveform between 0 and 1 that the oscillator slaves to
// * `res`: the frequency of resonance as a factor of the fundamental pitch.
//------------------------------------------------------------
// Author: Bart Brouns
// License: GPLv3
// CZ oscillators by Mike Moser-Booth:
// <https://forum.pdpatchrepo.info/topic/5992/casio-cz-oscillators>
// Ported from pd to Faust by Bart Brouns

CZresTrap(fund, res) = CZ.resTrap(fund, res);

CZ = environment {
  
    saw(fund, index) = sawChooseP(fund, index, 0);
    sawP(fund, index) = sawChooseP(fund, index, 1);
    sawChooseP(fund, index, p) =
      (((FUND(fund,align,p)*((.5-INDEX)/INDEX)),(-1*FUND(fund,align,p)+1)*((.5-INDEX)/(1-INDEX))):min+FUND(fund,align,p))*2*ma.PI:cos
    with {
      INDEX = (.5-(index*.5)):max(0.01):min(0.5);
      align = si.interpolate(index, 0.75, 0.5);
    };
    
    square(fund, index) = squareChooseP(fund, index, 0);
    squareP(fund, index) = squareChooseP(fund, index, 1);
    squareChooseP(fund, index, p) = (FUND(fund,align,p)>=0.5), (ma.decimal((FUND(fund,align,p)*2)+1)<:_-min(_,(-1*_+1)*((INDEX)/(1-INDEX)))) :+ *ma.PI:cos
    with {
      INDEX = (index:pow(0.25)):max(0):min(1);
      align = si.interpolate(INDEX, -0.25, 0);
    };

    pulse(fund, index) = pulseChooseP(fund, index, 0);
    pulseP(fund, index) = pulseChooseP(fund, index, 1);
    pulseChooseP(fund, index, p) = ((FUND(fund,align,p)-min(FUND(fund,align,p),((-1*FUND(fund,align,p)+1)*(INDEX/(1-INDEX)))))*2*ma.PI):cos
    with {
      INDEX = index:min(0.99):max(0);
      align = si.interpolate(index, -0.25, 0.0);
    };

    sinePulse(fund, index) = sinePulseChooseP(fund, index, 0);
    sinePulseP(fund, index) = sinePulseChooseP(fund, index, 1);
    sinePulseChooseP(fund, index, p) = (min(FUND(fund,align,p)*((0.5-INDEX)/INDEX),(-1*FUND(fund,align,p)+1)*((.5-INDEX)/(1-INDEX)))+FUND(fund,align,p))*4*ma.PI:cos
    with {
      INDEX = ((index*-0.49)+0.5);
      align = si.interpolate(index, -0.125, -0.25);
    };

    halfSine(fund, index) = halfSineChooseP(fund, index, 0);
    halfSineP(fund, index) = halfSineChooseP(fund, index, 1);
    halfSineChooseP(fund, index, p) = (select2(FUND(fund,align,p)<.5, .5*(FUND(fund,align,p)-.5)/INDEX+.5, FUND(fund,align,p)):min(1))*2*ma.PI:cos
    with {
      INDEX = (.5-(index*0.5)):min(.5):max(.01);
      align = si.interpolate(index:min(0.975), -0.25, -0.5);
    };
    
    FUND =
      case {
        (fund,align,0) => fund;
        (fund,align,1) => (fund+align) : ma.frac; // align phase with fund
      };
    resSaw(fund,res) = (((-1*(1-fund))*((cos((ma.decimal((max(1,res)*fund)+1))*2*ma.PI)*-.5)+.5))*2)+1;
    resTriangle(fund,res) = select2(fund<.5, 2-(fund*2), fund*2)*INDEX*2-1
    with {
	  INDEX = ((fund*(res:max(1)))+1:ma.decimal)*2*ma.PI:cos*.5+.5;
    };
    resTrap(fund, res) = (((1-fund)*2):min(1)*sin(ma.decimal(fund*(res:max(1)))*2*ma.PI));
};

//===============================PolyBLEP-Based Oscillators=================================

//----------`(os.)polyblep`----------
// PolyBLEP residual function - used for smoothing steps in the audio signal.
//
// #### Usage
//
// ```
// polyblep(Q, phase) : _
// ```
//
// Where:
//
// * `Q`: smoothing factor between 0 and 0.5. Determines how far from the ends of the phase interval the quadratic function is used.
// * `phase`: normalised phase (between 0 and 1)
//------------------------------------------------------------
// Author: Jacek Wieczorek
polyblep(Q, phase) = (0, L(phase / Q), R((phase - 1) / Q)) : select3(sel) 
with {
	sel = (phase < Q) + 2*(phase > 1 - Q);
	L(x) = 2*x - x*x - 1; // Used near the left end of the interval
	R(x) = 2*x + x*x + 1; // Used near the right end of the interval
};

//----------`(os.)polyblep_saw`----------
// Sawtooth oscillator with suppressed aliasing (using polyBLEP).
//
// #### Usage
//
// ```
// polyblep_saw(f) : _
// ```
//
// Where:
//
// * `f`: frequency in Hz 
//------------------------------------------------------------
// Author: Jacek Wieczorek
polyblep_saw(f) = naive - polyblep(Q , phase)
with {
	phase = phasor(1, f);
	naive = 2 * phase - 1;
	Q = f / ma.SR;
};

//----------`(os.)polyblep_square`----------
// Square wave oscillator with suppressed aliasing (using polyBLEP).
//
// #### Usage
//
// ```
// polyblep_square(f) : _
// ```
//
// Where:
//
// * `f`: frequency in Hz 
//------------------------------------------------------------
// Author: Jacek Wieczorek
polyblep_square(f) = naive - polyblep(Q, phase) + polyblep(Q, ma.modulo(phase + 0.5, 1))
with {
	phase = phasor(1, f);
	naive = 2 * (phase * 2 : int) - 1;
	Q = f / ma.SR;
};

//----------`(os.)polyblep_triangle`----------
// Triangle wave oscillator with suppressed aliasing (using polyBLEP).
//
// #### Usage
//
// ```
// polyblep_triangle(f) : _
// ```
//
// Where:
//
// * `f`: frequency in Hz 
//------------------------------------------------------------
// Author: Jacek Wieczorek
polyblep_triangle(f) = polyblep_square(f) : fi.pole(0.999) : *(4 * f / ma.SR);


//===============================Filter-Based Oscillators=================================

//-----------------`(os.)quadosc`--------------------
// Sinusoidal oscillator based on QuadOsc by Martin Vicanek.
//
// #### Usage
//
// ```
// quadosc(freq) : _
// ```
//
// where
//
// * `freq`: frequency in Hz
//
// #### Reference
// * <https://vicanek.de/articles/QuadOsc.pdf>
//------------------------------------------------------------
// Authors:
// Dario Sanfilippo <sanfilippo.dario@gmail.com>
// and Oleg Nesterov (jos ed.)
quadosc(f) = tick ~ (_,_) 
with {
    k1 = tan(f * ma.PI / ma.SR);
    k2 = 2 * k1 / (1 + k1 * k1);
    tick(u_0,v_0) = u_1,v_1 
    with {
        tmp = u_0 - k1 * v_0;
        v_1 = v_0 + k2 * tmp;
        u_1 = tmp - k1 * v_1 : select2(1',1);
    };
};

// end further contributions section
//########################################################################################