// Dedicated to Wendy Carlos who gave the world the first 
// glimpse of what can happen when anything is possible 
// in music.

#define kPI 3.1415926535897932384626433832795

global {

	srate 44100;
	krate  2100;
	interp 1;
	outchannels 2;

//	table tblSaw(buzz,2048, 7, 1, 0.5);
//	table tblBass(buzz,2048, 7, 1, -0.9);
	
	table tblSaw(harm,2048, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625);
	
// The same as buzz, but phase shifted so it doesn't cause
// quite so much trouble with interpolation. harm uses sum of sin;
// buzz uses sum of cos to construct a sharp spike of great harmonic 
// interest, but almost designed to break a table interpolator.
table tblBass(harm,2048, 
		0.15*1, 
		0.15*-0.9, 
		0.15*-0.9*-0.9, 
		0.15*-0.9*-0.9*-0.9, 
		0.15*-0.9*-0.9*-0.9*-0.9,
		0.15*-0.9*-0.9*-0.9*-0.9*-0.9,
		0.15*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9, 
		0.15*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9,
		0.15*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9,
		0.15*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9, 
		0.15*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9*-0.9
		);
		
	route (reverb_bus, AnalogSynth);
	route (reverb_bus, AnalogSynth2);
	route (reverb_bus, AnalogSynth3);
	route (reverb_bus, AnalogSynth4);
	send( Reverb; ; reverb_bus);
	route (output_bus, Reverb);

}
 
opcode combgain(xsig t, xsig rt) { 
  xsig temp; 

  temp = exp (log (10) * - 3 * t / rt); 
  return (temp); 
} 

aopcode xStereoReverb(asig inL, asig inR, ivar f0) {  
  return (inL*4, inR*4);
}


aopcode StereoReverb(asig inL, asig inR, ivar f0, ivar ...) {  
  ivar rt;
  ivar d2, distance, inv_distance, inv_d2;
  ksig first, ki; 
  asig ap1, ap2; 
  asig inMono;
  asig c[4]; 
  asig outL, outR;
  ivar t[4], revgain[4], i; 
  table timeconst (data, 5, 0.030000, 0.034300, 0.039300,
	                    0.045000, 0.000000); 
  oparray comb [4];

  rt = f0;
  
  inMono = (inL + inR)*0.5;
  distance = 0.600000; 
  inv_distance = 1.0/distance;
  d2 = distance * distance; 
  inv_d2 = 1.0/d2;
  i = 0; 
  while (i < 4) { 
    t [ i ] = tableread (timeconst, i); 
    revgain [ i ] = combgain(t [ i ], rt); 
    i = i + 1; 
    } 
  ap1 = allpass (inMono *inv_distance, 0.001700, 0.700000); 
  ap2 = allpass (ap1, 0.005000, 0.700000); 
  c[0]= comb (ap2, t[0], revgain[0]); 
  c[1]= comb (ap2, t[1], revgain[1]); 
  c[2]= comb (ap2, t[2], revgain[2]); 
  c[3]= comb (ap2, t[3], revgain[3]); 
  outL = (c[0]+ c[1]+ c[2]+ c[3]) *0.25; 
  outR = (c[0]- c[1]+ c[2]- c[3]) *0.25; 
  return (
      (inL *inv_d2 + outL) * (distance * 0.750000),
      (inR *inv_d2 + outR) * (distance * 0.750000)
      );
}

aopcode xResonantFilter(asig in, ivar idbResonance, ksig cut) {
  return (in);
}

aopcode ResonantFilter(asig in, ivar idbResonance, ksig cut)
{
	ivar Q;
	ivar alphaFactor;
	ksig omega;
	ksig sn, cs;
	ksig alpha;
	ksig b0, b1, b2;
	ksig a0, a1, a2;
	ksig b0inv;
	asig x0,x1;
	asig y0,y1;
	asig result;
	
	Q = pow(10, idbResonance*0.05);
	alphaFactor = 1/(2.0*Q);

	if (cut <= 0.0001) {
		// Roundoff conditions give the filter a DC drift 
		// as freq ----> 0, so guard against this.
		a0 = 0; a1 = 0; a2 = 0;
		b0 = 0; b1 = 0; b2 = 0;
	}else {
	
		omega = cut* (2.0*kPI/s_rate);
		sn = sin(omega);
		cs = cos(omega);
		
		// double alpha = sn/(2.0*Q);
		alpha = sn*alphaFactor;
	
		b0 = 1+alpha;
		b0inv = 1/(b0);
	
		
		a1 = (1-cs)*b0inv;
		a0 = (a1)*0.5;
		a2 = a0;
		
		b1 = (-2*cs)*b0inv;
		b2 = (1-alpha)*b0inv;
	}
	
	result = a0 *in+a1*x0 + a2*x1 - b1*y0 - b2*y1;

	x1 = x0;
	x0 = in;
	y1 = y0;
	y0 = result;
	return (result);
}


aopcode SawWave(xsig cps) {
	imports exports table tblSaw;
	asig oscil1;
	oscil1 = tableread(tblSaw, 2048*aphasor(cps));
	return (oscil1);
}
aopcode BassWave(xsig cps) {
	imports exports table tblBass;
	asig oscil1;
	oscil1 = tableread(tblBass, 2048*aphasor(cps));
	return (oscil1);
}


kopcode Lfo(ksig cps) {
	return ( sin (2*kPI * kphasor(cps)) );
}

instr Reverb() {
//   output(input*2);
  output (StereoReverb(input[0]*2, input[1]*2, 1.8, 0));
}



aopcode pan(asig in, ivar angle) 
{
	
	// equal-power panning with a delay
	ivar leftVol;
	ivar rightVol;
	ivar leftDelay;
	ivar rightDelay;
	ivar maxDelay;
	ivar phi;
	asig out[2];

	maxDelay = (0.1/33.0); // 0.2meters/(33 meters/second)
	phi = kPI*(angle+1);
	rightVol = sin(phi*0.5);
	leftVol = cos(phi*0.5);
	if (angle < 0) 
	{
		leftDelay = 1.0/s_rate;
		rightDelay = cos(-angle*(kPI*0.5))*maxDelay+1.0/s_rate;
	}  else {
		leftDelay = cos(angle*(kPI*0.5))*maxDelay+1.0/s_rate;
		rightDelay = 1.0/s_rate;
	}
	if (angle < 0) 
	{
		out[0] = delay(rightVol*in, rightDelay);
		out[1] =  leftVol*in;
	}  else {
		out[0] = rightVol*in;
		out[1] = delay(leftVol*in, leftDelay);
	}
	return (out);
}

kopcode ADSR(
	ivar attack, 
	ivar decayHalfTime,
	ivar sustainLevel, 
	ivar releaseHalfTime,
	ksig bReleased
	)
{
	ksig envPhase;
	ksig gain;
	ksig dGain;
	ksig pGain;
	ksig cycles;
	ksig retval;

	if (cycles <= 0 || bReleased)
	{
		if (envPhase == 0) { // Attack
			cycles = (attack*k_rate);
			gain = 0;
			dGain = 1/cycles;
			pGain = 1.0;
			envPhase = 1;
			if (cycles != 0) {
				return (0);
			}
		}
		cycles = 1000000000.0; // Just something very large.
		if (!bReleased) 
		{
			gain = 1.0;
			pGain = exp((log(0.5)/k_rate)/decayHalfTime);
			dGain = sustainLevel*(1-pGain);
		} else {
			pGain = exp((log(0.5)/k_rate)/decayHalfTime);
			dGain = 0;
		}
	}

	cycles = cycles-1;	
	gain = gain*pGain+dGain; // More expensive to compute, but less expensive than a branch.
	return (gain);

}

instr AnalogSynth(midiNote, velocity)
{
	asig out[2];
	asig aOscil;
	ksig kLfo;
	ksig kAdsr;
	ksig kFilterEnv;
	asig aResult;
	asig ringMod;
	ivar cps;
	imports ksig pitchbend; // event target.
	ksig cps_pb;

	kLfo = Lfo(5)*0.005 + 1;
	
	cps = cpsmidi(midiNote);
	cps_pb = cpsmidi(midiNote+pitchbend);
	aOscil = BassWave(cps_pb *kLfo);
	
	ringMod = BassWave(5.02*cps_pb )*10;
	aOscil = aOscil*ringMod+aOscil;
	
	kAdsr = ADSR(0.01, 0.1, 0.9, 0.01, released);
	kFilterEnv = kline(0, 0.05, cps*3, 0.1, cps*2, 0.5, cps, 100000, cps);
	
	
	
	if (released) {
		if (kAdsr > 0.00001) {
			extend(1/k_rate);
		}
	}
	aResult = ResonantFilter(velocity*aOscil *kAdsr *0.1, 6, kFilterEnv);
	out = (pan(aResult,0.7));
	output (out);
}

instr AnalogSynth2(midiNote, velocity)
{
	asig out[2];
	asig aOscil;
	ksig kLfo;
	ksig kAdsr;
	ksig kFilterEnv;
	asig aResult;
	ivar cps;

	kLfo = Lfo(5)*0.005 + 1;
	
	cps = cpsmidi(midiNote);
	aOscil = SawWave(cps*kLfo);
	kAdsr = velocity*ADSR(0.01, 0.1, 0.9, 0.01, released);
	kFilterEnv = kline(cps*5, 0.02, cps, 0.2, cps*0.5, 10000, cps*0.5);
//	kFilterEnv = kline(cps*5, 0.01, cps, 0.2, cps*0.5);
	
	if (released) {
		if (kAdsr > 0.00001) {
			extend(1/k_rate);
		}
	}
	aResult = ResonantFilter(aOscil *kAdsr *0.1, 6, kFilterEnv);
	out = (pan(aResult, (midiNote-72)*(1/3.5)));
	output (out);
}

instr AnalogSynth3(midiNote, velocity)
{
	asig out[2];
	asig aOscil;
	ksig kLfo;
	ksig kAdsr;
	ksig kFilterEnv;
	asig aResult;
	imports ksig pitchbend; // event target.
	ksig cps_pb;

	kLfo = Lfo(5)*0.005 + 1;
	
	cps_pb = cpsmidi(midiNote + pitchbend);
	
	aOscil = SawWave(cps_pb*kLfo);
	kAdsr = ADSR(0.15, 0.3, 1.0, 0.01, released);
	kFilterEnv = cps_pb*kline(0.001, 0.2, 1, 5, 5, 0.5, 1, 100000, 1);
	
	
	if (released) {
		if (kAdsr > 0.00001) {
			extend(1/k_rate);
		}
	}
	aResult = ResonantFilter(aOscil *kAdsr *(velocity*0.1/7.0), 12, kFilterEnv);
	out = (pan(aResult,-0.7));
	output (out);
}

instr AnalogSynth4(midiNote, velocity)
{
	asig out[2];
	asig aOscil;
	ksig kLfo;
	ksig kAdsr;
	ksig kFilterEnv;
	asig aResult;
	asig ringMod;
	ivar cps;
	imports ksig pitchbend; // event target.
	ksig cps_pb;

	kLfo = Lfo(5)*0.005 + 1;
	
	cps = cpsmidi(midiNote);
	cps_pb = cpsmidi(midiNote+pitchbend);
	aOscil = BassWave(cps_pb *kLfo);
	
	
	kAdsr = ADSR(0.01, 0.1, 0.9, 0.01, released);
	kFilterEnv = cps_pb*kline(0,  0.05, 5,  0.1,  2,  0.5, 3,  100, 1);
	
	
	
	if (released) {
		if (kAdsr > 0.00001) {
			extend(1/k_rate);
		}
	}
	aResult = ResonantFilter(aOscil *kAdsr *(velocity*0.1/7.0), 6, kFilterEnv);
	out = (pan(aResult,0.3));
	output (out);
}