#!/usr/bin/env perl

use Math::Complex;
use Math::Trig;

# proportion of total resistance for a log taper pot

sub logpot($) {
    my $x = shift;

    my $ym = 0.2;

    my $b = ( ( 1. / $ym ) - 1. );
    $b = $b * $b;

    my $a = 1. / ( $b - 1. );

    return $a * ( $b**$x ) - $a;
}

# proportion of resistance between two points on a log pot
sub logprop($$$) {
    my ( $v, $t1, $t2 ) = @_;

    $v  = logpot($v);
    $t1 = logpot($t1);
    $t2 = logpot($t2);

    return ( $v - $t1 ) / ( $t2 - $t1 );
}

# compute impedance of two elements in parallel
sub par($$) {
    my ( $a, $b ) = @_;

    return 1. / ( ( 1. / $a ) + ( 1. / $b ) );
}

# impedance of capacitor in nF at frequency f
sub cap($$) {
    my ( $c, $f ) = @_;

    my $z = $c * 1.e-9 * 2 * pi * $f;

    $z = -i / $z;

    return $z;
}

@fns  = ();
@frqs = ();

# calc gain at frequeny f and knob position k (0-10)
sub calc_gain($$) {
    my ( $f, $k ) = @_;

    # position of taps on pot.
    my $te = 0.4;    # terminal e
    my $td = 0.6;    # terminal d

    my $rac = 470000.;

    #terminal c is input
    #terminal b is output
    #terminal a is ground

    # loads on terminals d and e
    my $ld = 18000. + cap( 47., $f );    #18k+47nF
    my $le = 4700. + cap( 250., $f );    #4.7k+250nF


    #compute resistances of track
    my $rae = $rac * logpot($te);
    my $rad = $rac * logpot($td);

    my $red = $rad - $rae;

    my $rdc = $rac - $rad;

    my $tot = $rae + $red + $rdc;

    print "Track resistances: $rae + $red + $rdc  = $tot\n";

    # next compute the efective impedances of each leg

    my $erae = par( $rae,         $le );
    my $erad = par( $erae + $red, $ld );
    my $ered = $erad - $erae;
    my $erac = $erad + $rdc;
    my $erdc = $erac - $erad;

    my $etot = $erac;

    # input filter
    $etot += 150000. + cap( 15., $f );    # 150k +15nF
                                          #balance control
    $etot += par( 1100., cap( 120., $f ) );

    print "Effective total impedance $erae + $ered + $erdc = $erac\n";

    # compute the potentials for a 1V input

    my $pe  = $erae / $etot;
    my $pd  = $erad / $etot;
    my $pa  = $erac / $etot;
    my $ped = $pd - $pe;
    my $pda = $pa - $pd;

    my $v = $k / 10.;
    my $hf;

    if ( $v < $te ) {
        $hf = $pe * logprop( $v, 0, $te );
    }
    elsif ( $v < $td ) {
        $hf = $pe + ( $ped * logprop( $v, $te, $td ) );
    }
    else {
        $hf = $pd + ( $pda * logprop( $v, $td, 1. ) );
    }

    return abs($hf);
}

open P, "|gnuplot";
print P "set term png size 1280,1024\n";
print P "set output 'volume.png'\n";
print P "set logscale x\n";
print P "set logscale y\n";
print P "set yrange [:1.1]\n";
print P "set xlabel 'frequency/Hz'\n";
print P "set ylabel 'gain'\n";
print P "plot ";

for ( my $k = 1 ; $k < 10.5 ; ++$k ) {
    my $gmax = 0;
    for ( my $f = 20 ; $f < 20000.1 ; $f *= 1.1 ) {
        my $g = calc_gain( $f, $k );
        $gmax = $g if $g > $gmax;
    }
    my $fn = "volume-" . $k . ".dat";
    open F, ">", $fn;
    for ( my $f = 20 ; $f < 20000.1 ; $f *= 1.1 ) {
        my $g = calc_gain( $f, $k );
        print F $f, " ", $g, " ", $g / $gmax, "\n";
    }
    close F;

    print P " '" . $fn . "' using 1:2 with lines title '" . $k . "',";

}
print "\n";
close P;