#!/usr/bin/perl -w

# see antenna_simulation.readme for info.

use strict;

my $name = $0;
$name =~ s/\.pl$//;
$name =~ s|.*/||;
open(LOG, ">>$name.log");

sub pr {
    print LOG @_;
    print @_;
}

sub prf {
    printf LOG @_;
    printf @_;
}

sub eng($$$) {
    my $c = shift;
    my $sep = shift;
    my $unit = shift;
    my $str;

    if (0) {
    } elsif ($c >= 1e12) {
	$str = sprintf "%.2f%sT%s", $c*1e-12, $sep, $unit;
    } elsif ($c >= 1e9) {
	$str = sprintf "%.2f%sG%s", $c*1e-9, $sep, $unit;
    } elsif ($c >= 1e6) {
	$str = sprintf "%.2f%sM%s", $c*1e-6, $sep, $unit;
    } elsif ($c >= 1e3) {
	$str = sprintf "%.2f%sk%s", $c*1e-3, $sep, $unit;
    } elsif ($c >= 1) {
	$str = sprintf "%.2f%s%s", $c, $sep, $unit;
    } elsif ($c >= 1e-3) {
	$str = sprintf "%.2f%sm%s", $c*1000, $sep, $unit;
    } elsif ($c >= 1e-6) {
	$str = sprintf "%.2f%su%s", $c*1e6, $sep, $unit;
    } elsif ($c >= 1e-9) {
	$str = sprintf "%.2f%sn%s", $c*1e9, $sep, $unit;
    } elsif ($c >= 1e-12) {
	$str = sprintf "%.2f%sp%s", $c*1e12, $sep, $unit;
    } else {
	$str = $c . $sep . $unit;
    }
    $str;
}

my $n = 2;
if (defined($ARGV[0])) { $n = $ARGV[0]; }

my $pi = 3.14159265358979323844;
my $f = 13.56e6;
my $l = 440e-9 / 4 * $n * $n;
my $w = 2*$pi*$f;
my $w2 = $w*$w;
my $ctot = 1/($w2*$l);

my $cpar = 10e-12; # midpoint of a 2-18pF variable capacitor
if (defined($ARGV[1])) { $cpar = $ARGV[1] * 1e-12; }
my $rpar = 2*$pi*$f*$l*20;
if (defined($ARGV[2])) { $rpar = $ARGV[2]; }

my $c = $ctot - 10e-12;
my $c2 = $c * sqrt($rpar / 50);
if (defined($ARGV[3])) { $c2 = $ARGV[3] * 1e-12; }

my $c1 = 1/(1/$c - 1/$c2);
my $c1a = 2*$c1;
if (defined($ARGV[4])) { $c1a = $ARGV[4] * 1e-12; }

pr "Cmdline N Cpar Rpar C2 C1a:", join(" ", @ARGV), "\n";

open(FH, ">Todo/$name.ckt") || die("cannot write to $name.ckt");
print FH "$name

V1   1 0 generator(1) ac 1
Ro   1 2 50
C1a  2 4 $c1a
C1b  0 3 $c1a
C2   2 0 $c2
Cpar 4 3 $cpar
L    4 3 $l
Rpar 4 3 $rpar

.option nopage
.print ac I(L) V(2)
.AC 12.5Meg 14.5Meg 10k > Todo/$name.dat
.end
";
close FH;

pr "N:    $n
L:    $l
Cpar: $cpar
Rpar: $rpar
Ctot: $ctot
C2:   $c2
C1a:  $c1a
";

`gnucap -b Todo/$name.ckt >/dev/null`;

open(FH, "Todo/$name.dat") || die("cannot open Todo/$name.dat");
my $save;
my $max = 0;
my @data;
my $ix = 0;
while (<FH>) {
    my $line = $_;

    next if (m/^#/);
    tr/ \t/ /s;
    s/^ //;
    my @fld = split;
    push @data, [ @fld ];
    if ($max < $fld[1]) {
	$save = $ix;
	$max = $fld[1];
    }
    $ix++;
}

my $lvl = $max * sqrt(0.5);
my $f1 = $data[0][0];
my $f2 = $data[$#data][0];

sub prt($$) {
    my $prefix = shift;
    my $ix = shift;
    my $r = $data[$ix];
    printf LOG "%s%.3f %s %s\n", $prefix, $$r[0]/1e6, $$r[1], $$r[2];
    printf "%s%.3f %s %s\n", $prefix, $$r[0]/1e6, $$r[1], $$r[2];
}

my $i1;
my $i2;
for ($ix = 0; $ix < $save; $ix++) {
    if ($data[$ix][1] >= $lvl ) {
	$f1 = $data[$ix][0];
	prt("f1:  ", $ix);
	$i1 = $ix;
	last;
    }
}
prt("max: ", $save);

#print $save, " ", $#data, "\n";
for ($ix = $save + 1; $ix < @data; $ix++) {
    if ($data[$ix][1] <= $lvl ) {
	$f2 = $data[$ix][0];
	prt("f2:  ", $ix);
	$i2 = $ix;
	last;
    }
}

my $q = $data[$save][0] / ($f2 - $f1);
prf "Q = %.1f\n", $q;

$l    = eng($l,    " ", "H");
$cpar = eng($cpar, " ", "F");
$rpar = eng($rpar, " ", "Ohm");
$ctot = eng($ctot, " ", "F");
$c2   = eng($c2,   " ", "F");
$c1a  = eng($c1a,  " ", "F");
$q = sprintf "%.1f", $q;

open(FH, ">Todo/$name.plot") || die("cannot write to $name.plot");
print FH "
set grid xtics ytics
set tics out

set xlabel \"MHz\"
set ylabel \"Current through loop [mA]\"
set y2label \"Input voltage [mV]\"

set yrange [0:100]
set ytics 10

set y2range [0:1000]
set y2tics 100

set label \"N\"     at graph 0.03,0.95
set label \"L\"     at graph 0.03,0.92
set label \"Cpar\"  at graph 0.03,0.89
set label \"Rpar\"  at graph 0.03,0.86
set label \"Ctot\"  at graph 0.03,0.83
set label \"C2\"    at graph 0.03,0.80
set label \"C1a\"   at graph 0.03,0.77
set label \"Q\"     at graph 0.03,0.74

set label \"$n\"    at graph 0.08,0.95
set label \"$l\"    at graph 0.08,0.92
set label \"$cpar\" at graph 0.08,0.89
set label \"$rpar\" at graph 0.08,0.86
set label \"$ctot\" at graph 0.08,0.83
set label \"$c2\"   at graph 0.08,0.80
set label \"$c1a\"  at graph 0.08,0.77
set label \"$q\"    at graph 0.08,0.74

";

my $fm = $data[$save][0]*1e-6;
my $fm1 = $data[$save][1]*1000 + 6;
my $fm2 = $data[$save][1]*1000 + 1;

print FH "set arrow from $fm,$fm1 to $fm,$fm2\n";

if (defined($f1)) {
    $f1 = $data[$i1][0]*1e-6;
    my $f11 = $data[$i1][1]*1000 + 6;
    my $f12 = $data[$i1][1]*1000 + 1;
    print FH "set arrow from $f1,$f11 to $f1,$f12\n";
}
if (defined($f2)) {
    $f2 = $data[$i2][0]*1e-6;
    my $f21 = $data[$i2][1]*1000 + 6;
    my $f22 = $data[$i2][1]*1000 + 1;
    print FH "set arrow from $f2,$f21 to $f2,$f22\n";
}

print FH "
set terminal postscript color enhanced 10 # eps font \"Courier,12\"
set output \"Todo/$name.ps\"

plot \\
 \"Todo/$name.dat\" using (\$1*1e-6):(\$2*1000) with lines lt 1 lc rgb \"black\" title \"I(L)\",\\
 \"Todo/$name.dat\" using (\$1*1e-6):(\$3*1000) with lines lt 5 lc rgb \"black\" title \"V(2)\" axes x1y2
";
close FH;

`gnuplot Todo/$name.plot`;

__END__

sub cadd(@) {
    my $s = 0;
    for my $ci (@_) {
	$s += 1/$ci;
    }
    1/$s;
}
