#!/usr/bin/perl -w

use strict;
use Opt;
use Prefix;
use Eseries;

#use Math::Trig ':pi';
use Math::Complex qw(:trig :pi);
use Data::Dumper;

#############

my $sqrt2 = sqrt(2);
#my $pi = 2 * atan2(1,0);
my $log10 = log(10);
my %disp = ();

my $Usage;
sub Usage(;$) {
    my $str = shift;

    if ($str) {
	print $str;
	print "\n";
    }
    my $prg = $0;
    $prg =~ s|.*/||;

    print
"Usage:
\t$prg <function> N=.. key=val ...
Available functions:
";
    for my $k (sort keys %disp) {
	my $r = $disp{$k};
	my ($func, $doc) = @$r;
	printf "\t%-15s %s\n", "$k:", $doc;
    }

    if ($str) {
	exit 1;
    }
    exit 0;
}

sub options_defaults($) {
    my $data = shift;
    %$data = (
	N  => 1,
    );
}
sub options_check($) {
    my $data = shift;

    my $N = $$data{N};

    Opt::err_eval( undef, ">  0", qw/ N /);
    if ($N == int($N)) {
	# ok
    } else {
	Usage("N ($N) must be a positive integer");
    }
}
sub act_filter_init($;$) {
    my $data = shift;
    $Usage = shift;

    options_defaults($data);
    @ARGV = Opt::opt($data, 0, @ARGV);
    Prefix::decode_hash($data);
    options_check($data);
    $data;
}

#############

sub exp10($) {
    # 10^x
    my $x = shift;
    my $y = exp($x*$log10);
    $y
}
sub dbToKvot($) {
    my $x = shift;
    my $y = exp10($x/20);
    $y
}

sub f_rel($) {
    my $data = shift;
    my ($a, $b) = @$data;
    my $f_rel = 1;

    if ($b != 0) {
	my $c = 2*$b - $a*$a;
	$f_rel = sqrt( $c + sqrt( $c * $c + 4*$b*$b ) ) / ( $sqrt2 * $b);
    }

    $f_rel;
}

sub Q($) {
    my $data = shift;
    my ($a, $b) = @$data;
    my $Q = sqrt($b)/$a;

    $Q;
}

sub row(@) {
    my $type = shift;
    my @data = @_;


    if ($type eq "ab") {
	printf(" a        b         f/fc     fc/f     Q        1/Q       1/sqrt(b)\n");
	for my $r (@data) {
	    my ($a, $b) = @$r;
	    my $f_rel = f_rel($r);
	    my $Q = Q($r);
	    my $Qi = $Q == 0 ? "9" : 1/$Q;
	    if ($b == 0) {
		printf "%.6f %.6f  %.6f %.6f %.6f %.6f\n", $a, $b, $f_rel, 1/$f_rel, $Q, $Qi;
	    } else {
		printf "%.6f %.6f  %.6f %.6f %.6f %.6f  %.6f\n", $a, $b, $f_rel, 1/$f_rel, $Q, $Qi, 1/sqrt($b);
	    }
	}
    } elsif ($type eq "c") {
	for my $r (@data) {
	    printf " %.6f", $r;
	}
	print "\n";
    }
}

sub poly_add($$) {
    my $a = shift;
    my $b = shift;

    my $la = @$a;
    my $lb = @$b;

    my $ll = $la < $lb ? $lb : $la ;
    my @c = (0) x $ll;
    for (my $ix = 0; $ix < $ll; $ix++) {
	my $ae = $$a[$ix] // 0;
	my $be = $$b[$ix] // 0;
	$c[$ix] = $ae + $be;
    }
    [ @c ];
}
sub poly_mul($$) {
    my $a = shift;
    my $b = shift;

    my $la = @$a;
    my $lb = @$b;

    my $ll = $la + $lb;
    my @c = (0) x $ll;
    for (my $xa = 0; $xa < $la; $xa++) {
	for (my $xb = 0; $xb < $lb; $xb++) {
	    my $ea = $$a[$xa] // 0;
	    my $eb = $$b[$xb] // 0;
	    $c[$xa + $xb] += $ea * $eb;
	}
    }
    [ @c ];
}

##############

sub kriticher($) {
    my $data = shift;
    my $N = $$data{N};
    if ($N < 1 || int($N) != $N) {
	warn("N = $N should be a positive integar");
	return;
    }

    my $alfa = sqrt( exp( (1/$N) * log(2) ) - 1 );

    my $b = $alfa * $alfa;
    my $a = 2* $alfa;
    my @data = ("ab");
    my $cnt = int($N/2);
    my $rst = $N - 2*$cnt;

    if ($rst) {
	push @data, [ $alfa, 0 ];
    }
    push @data, ([ $a, $b ])x$cnt;

    #print "$alfa $a $b\n";

    @data;
}
$disp{kritich} = [ \&kriticher, "critcally damped filter coeff. at order N" ];

sub bessel_polynomial($) {
    my $data = shift;
    my $N = $$data{N};
    my $c = 1;

    # coefficients for non factorized polynom
    # A, with c[0] == 1
    my @c = ("c", 1);
    for (my $ix = 1; $ix <= $N; $ix++) {
	$c *= 2*($N - $ix + 1)/($ix * ( 2*$N - $ix + 1));
	#printf "%.6f\n", $c;
	push @c, $c;
    }

    # B, with c[N] = 1
    my $k = $c[$#c];
    for (my $ix = 1; $ix < @c; $ix++) {
	$c[$ix] /= $k;
	#print "$c[$ix]\n";
    }
    #print join(" ", @c), "\n";

    @c;
}
$disp{bessel_polynomial} = [ \&bessel_polynomial, "Bessel polynomial" ];

sub bessel($) {
    my $data = shift;
    my $N = $$data{N};

    if ($N < 1 || int($N) != $N) {
	warn("N = $N should be a positive integar");
	return;
    }
    my $b = (-1 + sqrt(5))/2;
    my $a = sqrt(3*$b);

    #printf "%.6f %.6f %.6f\n", $a/$b, 1/$b, 1/sqrt($b);

    # from Halbleiterschaltungstechnik, U. Tietze, Ch. Schenk, E. Gamm, Springer
    my @data = (
	[ [ 1, 0 ], ],
	[ [ $a, $b ] ],
	[ [ 0.7560, 0.0000 ], [ 0.9996, 0.4472 ] ],
	[ [ 1.3397, 0.4889 ], [ 0.7743, 0.3890 ] ],
	[ [ 0.6656, 0.0000 ], [ 1.1402, 0.4128 ], [ 0.6216, 0.3245 ] ],
	[ [ 1.2217, 0.3887 ], [ 0.9686, 0.3505 ], [ 0.5131, 0.2856 ] ],
	[ [ 0.5937, 0.0000 ], [ 1.0944, 0.3395 ], [ 0.8304, 0.3011 ], [ 0.4332, 0.2381 ] ],
	[ [ 1.1112, 0.3162 ], [ 0.9754, 0.2979 ], [ 0.7202, 0.2621 ], [ 0.3728, 0.2087 ] ],
	[ [ 0.5386, 0.0000 ], [ 1.0244, 0.2834 ], [ 0.8710, 0.2636 ], [ 0.6320, 0.2311 ], [ 0.3257, 0.1854 ] ],
	[ [ 1.0215, 0.2650 ], [ 0.9393, 0.2549 ], [ 0.7815, 0.2351 ], [ 0.5604, 0.2059 ], [ 0.2883, 0.1665 ] ],
    );

    my @res = @{$data[$N-1]};
    unshift @res, "ab";
    @res;
}
$disp{bessel} = [ \&bessel, "Bessel approx. tabled values" ];

sub butterworth($) {
    my $data = shift;
    my $N = $$data{N};
    if ($N < 1 || int($N) != $N) {
	warn("N = $N should be a positive integar");
	return;
    }
    my @data = ("ab");

    if (2*int($N/2) == $N) {
	# even N
	for (my $ix = 1; $ix <= $N/2; $ix++) {
	    my $a = 2 * cos( (2*$ix - 1) * pi / (2*$N) );
	    my $b = 1;
	    push @data, [ $a, $b ];
	}
    } else {
	# odd N
	push @data, [ 1, 0 ];
	for (my $ix = 2; $ix <= ($N+1)/2; $ix++) {
	    my $a = 2 * cos( ($ix - 1) * pi / $N );
	    my $b = 1;
	    push @data, [ $a, $b ];
	}
    }
    @data;
}
$disp{butterworth} = [ \&butterworth, "Butterworth approx. coeff." ];

sub tjebytjov($) {
    my $data = shift;
    my $N = $$data{N};

    if ($N < 1 || int($N) != $N) {
	warn("N = $N should be a positive integar");
	return;
    }
    my $db = $$data{db} // 3;
    my $kvot = dbToKvot($db);
    my $epsilon = sqrt($kvot*$kvot -1);
    my $ny = asinh(1/$epsilon)/$N;
    my $aa = cosh( acosh(1/$epsilon)/$N ); # omega at 3dB limit
    my $a2 = $aa*$aa;

    my @data = ("ab");

    my $b1 = cosh($ny);
    my $a1 = sinh($ny);
    if (2*int($N/2) == $N) {
	# even N
	for (my $ix = 1; $ix <= $N/2; $ix++) {
	    my $b2 = cos((2*$ix - 1)*pi/(2*$N) );
	    my $b = 1/($b1*$b1 - $b2*$b2);
	    my $a = 2 * $b * $a1 * cos( (2*$ix - 1) * pi / (2*$N) );
	    push @data, [ $a * $aa, $b * $a2 ];
	}
    } else {
	# odd N
	push @data, [ 1/sinh($ny), 0 ];
	for (my $ix = 2; $ix <= ($N+1)/2; $ix++) {
	    my $b2 = cos(($ix - 1)*pi/($N) );
	    my $b = 1/($b1*$b1 - $b2*$b2);
	    my $a = 2 * $b * $a1 * cos( ($ix - 1) * pi / ($N) );
	    push @data, [ $a * $aa, $b * $a2 ];
	}
    }
    @data;
    ## TODO: the coefficients doesn't match table
    ## I guess $aa is wrong.
}
$disp{tjebytjov} = [ \&tjebytjov, "Tjebytjov approx. coeff." ];

sub tjebytjov_polynomial($) {
    my $data = shift;
    my $N = $$data{N};

    if ($N < 0 || int($N) != $N) {
	warn("N = $N should be a positive integar");
	return;
    }

    my @data = ("c");

    my @t = (
	[ 1 ],
	[ 0, 1 ],
    );
    for (my $ix = 2; $ix <= $N; $ix++) {
	my $p = poly_mul([ 0, 2 ], $t[$ix-1]);
	my $q = poly_mul([ -1 ], $t[$ix-2]);
	push @t, poly_add($p, $q);
    }
    my @c = @{$t[$N]};
    for (my $ix = $#c; $ix >= 0; $ix-- ) {
	last if ($c[$ix]);
	pop @c;
    }

    push @data, @c;

    @data;
}
$disp{tjebytjov_polynomial} = [ \&tjebytjov_polynomial, "Tjebytjov polynomial" ];

##############

sub mfb($) {
    my $data = shift;
    my $N = $$data{N} // 2;
    my $stage = $$data{stage} // 1;

    my $approx = $$data{approx};
    my $ref = $disp{$approx};
    if (!defined($ref)) {
	Usage("ERR: Unknown approximation: $approx\n");
    }
    my ($func,$doc) = @{$ref};

    my @arr = &$func($data);

    my ($a, $b) = @{$arr[$stage]};

    my $A0 = $$data{A0} // -1;
    my $fg = $$data{fg};
    my $C1 = $$data{C1};
    my $C2 = $$data{C2};
    my $Clim = undef;
    my $E  = $$data{E};
    my $R2 = $$data{R2};

    if (defined($C1) && defined($C2)) {
	$Clim = undef; # ???
    } elsif (defined($C1)) {
	$Clim = 4*$b * (1 - $A0) * $C1 / ($a*$a);
	my ($Cp, $Cc, $Cn)  = Eseries::nearest($Clim, @Eseries::E6);
	$C2 = $Cn;
    } elsif (defined($C2)) {
	$Clim = $C2 * $a * $a / (4*$b * (1 - $A0));
	my ($Cp, $Cc, $Cn)  = Eseries::nearest($Clim, @Eseries::E6);
	$C1 = $Cp;
    } else {
	Usage("Need C1 and/or C2\n");
    }

    if (defined($fg)) {
	$R2 = ( $a * $C2 - sqrt( $a*$a*$C2*$C2 - 4*$C1*$C2*$b*(1-$A0)) ) / (4*pi*$fg* $C1 * $C2);
    } elsif (defined($R2)) {
	$fg = ( $a * $C2 - sqrt( $a*$a*$C2*$C2 - 4*$C1*$C2*$b*(1-$A0)) ) / (4*pi*$R2* $C1 * $C2);
    } else {
	Usage("Need fg or R2\n");
    }

    $$data{A0} = $A0;
    $$data{C1} = $C1;
    $$data{fg} = $fg;
    $$data{Clim} = $Clim;
    $$data{C2} = $C2;
    $$data{R1} = - $R2 / $A0;
    $$data{R2} = $R2;
    $$data{R3} = $b / (4*pi*pi*$fg*$fg* $C1 * $C2 * $R2);

    if (defined($E)) {
	my @E = Eseries::find($E);
	my ($prev, $next);
	($prev, $$data{R1}, $next)  = Eseries::nearest($$data{R1}, @E);
	($prev, $$data{R2}, $next)  = Eseries::nearest($$data{R2}, @E);
	($prev, $$data{R3}, $next)  = Eseries::nearest($$data{R3}, @E);
    }

    #Prefix::println(1,[$data], " ",qw/ C1 F Clim F C2 F /);
    Prefix::println(1,[$data], " ",qw/ C1 F C2 F fg Hz/);
    Prefix::println(1,[$data], " ",qw/ R1 Ohm R2 Ohm R3 Ohm/);

}
$disp{mfb} = [ \&mfb, "multiple feedback circuit, use approx= A0= fg= C1= C2= stage=" ];

sub main() {
    my %data;

    my $data = act_filter_init(\%data, \&Usage);
    %data = %$data;

    my $function = shift @ARGV;
    if (defined($function)) {
    } else {
	Usage("");
    }

    my $ref = $disp{$function};
    if (!defined($ref)) {
	Usage("ERR: Unknown function: <$function>\n");
    }
    my ($func,$doc) = @{$ref};
    my @data = &$func(\%data);
    row(@data);
}

main;