Phaser modulation

It appears my math for the modulation function on the all-pass filter is not the best model for a JFET phaser. I spent a good couple hours with my guitar trying different settings on this phaser. There are certainly some nice sounds available, but it simply did not behave as expected.

I sat down and computed the bilinear transform myself, to better understand the sources I was reading, as well as doing my own analysis on small-signal resistance characteristic of a JFET. After 3 pages of algebra, I came to this conclusion:

R = 250/(1.0114 - sqrt(Vm));

This constrains the range of resistance to 250 to 22k-ohms. The constants 1.0114 and 250 include effects of parameters derived from 2N5458 JFET datasheet, also taking into account the parallel 22k resistor found in the MXR phase90 circuit.

You may notice that when put into the 1/2*pi*R*C equation for center frequency, we have an approximately linear frequency sweep up until the last bit. That explains some of the nature of the analog phaser. To give it a little bit of an exponential feel, we could do this:
R = 250/(lfo^2 + .0114);
or
R = 250/(0.5*(lfo + lfo^2) + .0114); ...to naturalize it a little bit. These are less accurate approximations, but are also less CPU intensive calculations than sqrt().

The final equation for "gain" in the difference equation (gl and gr in rakarrack code) is this:

g = [2*fs*C - (1.0114 - sqrt(Vm))/250 ] / [2*fs*C - (1.0114 - sqrt(Vm))/250 ]

Where fs is sampling frequency and C is value of capacitor (.05 uF).

The difference equation is this:
y[n] = g*(x[n] + y[n-1]) - x[n-1] where g is as noted above

Now a real JFET doesn't have a perfectly linear resistance characteristic, so we can introduce an input-dependent term that changes the "resistance" with signal amplitude:

Vm = lfo * (1-d*Vin^2) where constant "d" is amount of distortion to apply.

Vin is the input signal (smps[i]).

I think a realistic value for "d" is about 0.04 in this application. This can be a parameter that is possible to adjust from 0 to an insane number like .5 or .8 where there is some definite intermodulation available. This is not a rigorous account for the distortion added by the JFET, but I'm hoping it's a good enough approximation for acceptable sounding results while keeping CPU requirement at a minimum. The point is to add a small amount of harmonics at each phase stage to accentuate the "swoosh" heard from an analog phaser.

I see need to employ a look-up table for the LFO, so the phaser doesn't need to be interpolating it all the time.