function [R, V2mag, V2ang, K, freq_rng] = VCVSfilt3(C, fc, freq_rng, pass, type, order, plot_flg)
% rev 3: 12.18.03 correct fn on HP : thx to Ben Sprung
% 11.6.01 rev 2
% [R, V2mag, V2ang, K, freq_rng] = VCVSfilt2(C, fc, freq_rng, pass, type, order, plot_flg);
% en123 filter function, second order
% H & H chpt 5, 10.4.01
% C is capacitance 0.3e-6, for example
% R is R1 = R2 from H&H design, not the gain-setting resistor R = 10K
% fc is center frequency Hz, pass: 1 is LP 2 is HP 3 is BP
% if freq_rng = 0 then an internal range is created...
% type is butterworth / chebyshev 2dB / bessel == 1 2 3
% if plot_flg == 0 1 2 ==> no plot / mag / only mag+ang
% below: example of butterworth bandpass from 1000 to 10000 freq/sec, over range
% >> freq_rnge = 100:100:200000;
% >> figure, [RL, V2magL, V2ang, KL, freq_rangeL] = VCVSfilt3(1e-6, 10000, freq_rnge, 1, 1, 2, 1);
% >> figure, [RH, V2magH, V2ang, KH, freq_rangeH] = VCVSfilt3(1e-6, 5000, freq_rnge, 2, 1, 2, 1);
% >> figure, loglog(freq_rnge, V2magL .* V2magH)
% use + for bandreject

omeg_c = fc*2*pi;

if type == 2
   if pass == 1
      fc = .907*fc; % butterworth stays at fc alone
   elseif pass == 2 % HP
      fc = (1/.907)*fc;
   end
elseif type == 3
   if pass == 1
   fc = 1.272*fc;
   elseif pass == 2
      fc = (1/1.272)*fc;
   end
end

R = 1/(2*pi*fc*C);

if freq_rng == 0
   omeg_rng = round(omeg_c/20):ceil(omeg_c/50):round(omeg_c*20);
else
   omeg_rng = 2*pi*freq_rng;
end

or_sze = size(omeg_rng, 2);
Gvec = (1/R)*ones(1,or_sze);

if pass == 1
   Y1 = Gvec;
   Y2 = j*omeg_rng*C;
   Y3 = Gvec;
   Y4 = j*omeg_rng*C;
elseif pass == 2
   Y1 = j*omeg_rng*C;
   Y2 = Gvec;
   Y3 = j*omeg_rng*C;
   Y4 = Gvec;
elseif pass == 3
   Y1 = Gvec;
   Y2 = Gvec;
   Y3 = j*omeg_rng*C;
   Y4 = Gvec + j*omeg_rng*C;
end

if type == 1
   K = 1.586;
elseif type == 2
   K = 2.114;
elseif type == 3
   K = 1.268;
end

% forming the node admittance matrix naY by inspection:
% assume inward current to a node is negative.
% because of the voltage dependent voltage source
% at the op amp output, a -K*Y2 term is needed in
% naY position 1,2, making the off diagonal asymmetric...

for cnt = 1:size(omeg_rng,2)
    naY = [ Y1(cnt)+Y2(cnt)+Y3(cnt) -K*Y2(cnt)-Y3(cnt) ; -Y3(cnt) Y3(cnt)+Y4(cnt) ];
   % node admittance matrix
   Is = [1*Y1(cnt); 0]; % source vector

   V = inv(naY)*Is; % V = Y inverse * Is, the answer vector, [V1; V2]
   V2mag(cnt) = K*abs(V(2,:)); % output is K times larger than V2, or V+
   V2ang(cnt) = angle(V(2,:));
end

if freq_rng == 0 freq_rng = omeg_rng/(2*pi); end

if plot_flg == 1 | plot_flg == 2
   loglog(freq_rng, V2mag, '.')
   title('loglog spectrum of filter magnitude vs freq in Hz')
end
if plot_flg == 2
   figure
   semilogx(freq_rng, V2ang, 'r.')
   title('loglog spectrum of filter angle vs freq in Hz')
end