#!/opt/ActiveTcl-8.6/bin/tclsh
# Hey Emacs, use -*- Tcl -*- mode

set scriptname [file rootname $argv0]


# ---------------------- Command line parsing -------------------------
package require cmdline
set usage "usage: [file tail $argv0] \[options] filename"
set options {
}

# Bring in complex number functions
try {
    set complexnumbers_version [package require math::complexnumbers]
    puts "Loaded math::complexnumbers version $complexnumbers_version"
    # Always use the complex number versions of operands
    namespace import math::complexnumbers::*
} trap {TCL PACKAGE UNFOUND} {message optdict} {
    # message will be the human-readable error.  optdict will be the
    # options dictionary present when the command was attempted.
    puts "$message"
    puts [dict get $optdict -errorcode]
    exit
}

try {
    array set params [::cmdline::getoptions argv $options $usage]
} trap {CMDLINE USAGE} {msg o} {
    # Trap the usage signal, print the message, and exit the application.
    # Note: Other errors are not caught and passed through to higher levels!
    puts $msg
    exit 1
}


# ----------------------- Gnuplot settings ----------------------------

# The wxt terminal can keep windows alive after the gnuplot process
# exits.  This allows calling multiple persistent windows which allow
# zooming and annotation.
set gnuplot_terminal wxt


set pi 3.1415926535897932385


proc intrange {start points} {
    # Return a list of integers starting with start  with length points
    set count 1
    set integer $start
    set intlist [list]
    while {$count <= $points} {
	lappend intlist $integer
	incr integer
    }
    return $intlist
}

proc logrange {start perdecade end} {
    # Return a logarithmic points-long list of numbers starting at
    # start, having perdecade points per decade, and ending at end
    set alpha [expr 10 ** (double(1)/$perdecade)]
    set listnum $start
    set loglist [list]
    while {$listnum < $end} {
	lappend loglist $listnum
	set listnum [expr $listnum * $alpha]
    }
    lappend loglist $end
    return $loglist
}

proc magnum {complex} {
    # Returns the magnitude of the complex number
    set num [real [sqrt [* $complex [conj $complex]]]]
    return $num
}

proc zcoil {frequency re fc qmc qec} {
    # Return the (complex) coil impedance at a given frequency
    #
    # Arguments:
    #
    # frequency -- frequency in Hz
    # re -- DC resistance
    # fc -- resonant frequency (Hz)
    # qmc -- Mechanical Q
    # qec -- Electrical Q
    global pi
    set wc [expr 2 * $pi * $fc]
    # Separate things into terms to make debugging easier
    set wqt [expr $wc/$qmc + $wc/$qec]
    set wqm [expr $wc/$qmc]
    set s_cpx [complex 0 [expr 2 * $pi * $frequency]]
    set s2_cpx [* $s_cpx $s_cpx]
    set sqt_cpx [* $s_cpx [complex $wqt 0]]
    set w2_cpx [complex [expr $wc ** 2] 0]
    set sqm_cpx [* $s_cpx [complex $wqm 0]]
    set num12_cpx [+ $s2_cpx $sqt_cpx]
    set num_cpx [+ $num12_cpx $w2_cpx]
    set den12_cpx [+ $s2_cpx $sqm_cpx]
    set den_cpx [+ $den12_cpx $w2_cpx]
    set normalized_cpx [/ $num_cpx $den_cpx]
    set zcoil_cpx [* [complex $re 0] $normalized_cpx]
    return $zcoil_cpx
}

proc write_gnuplot_data {filename freq_list mag_list} {
    # Write a data file that can be plotted by gnuplot
    #
    # Arguments:
    #   filename -- Name of the output data file
    #   freq_list -- List of frequencies
    #   mag_list -- List of impedance magnitudes
    global params
    if { [catch {open $filename w} fp] } {
	puts "Could not open $filename for writing"
	return
    }
    foreach freq $freq_list mag $mag_list {
	puts $fp "$freq $mag"
    }
    close $fp
}

proc write_gnuplot_script {filename datafile} {
    # Write a gnuplot script to plot data in datafile
    #
    # Arguments:
    #   filename -- Name of the gnuplot script to be written
    #   datafile -- Name of the file containing data to be plotted
    global gnuplot_terminal
    global scriptname
    global speaker_dict
    if { [catch {open $filename w} fp] } {
	puts "Could not open $filename for writing"
	return
    }
    # Set the fractional amount of padding to require at the top
    # of the y-axis.  This makes more room for the plot legend.  A
    # value of 1 will give no padding, and a value of 2 will make
    # the padding twice the maximum current value plotted.
    set plot_padding 1.01
    # Annotation x postion in graph (0 --> 1) coordinates
    set anno_x 0.6
    # Annotation y position in graph (0 --> 1) coordinates
    set anno_y 0.5
    # Spacing between annotation lines in graph (0 --> 1) coordinates
    set anno_incr 0.1
    puts $fp "reset"
    puts $fp "set terminal $gnuplot_terminal size 800,600"
    puts $fp "set format y '%0.0s %c'"
    puts $fp "set format x '%0.0s %c'"
    puts $fp "set xlabel 'Frequency (Hz)'"
    puts $fp "set ylabel '|Z| (ohms)'"
    puts $fp "set nokey"
    puts $fp "set title '$scriptname'"
    # Start annotations
    puts $fp "set label 'fc = [dict get $speaker_dict fc] Hz' at graph $anno_x,$anno_y"
    set anno_y [expr $anno_y - $anno_incr]
    puts $fp "set label 'Qmc = [dict get $speaker_dict qmc]' at graph $anno_x,$anno_y"
    set anno_y [expr $anno_y - $anno_incr]
    puts $fp "set label 'Qec = [dict get $speaker_dict qec]' at graph $anno_x,$anno_y"
    set anno_y [expr $anno_y - $anno_incr]
    puts $fp "set label 'Re = [dict get $speaker_dict re] Ohms' at graph $anno_x,$anno_y"
    set anno_y [expr $anno_y - $anno_incr]
    set plot_instruction "using 1:2 with lines"
    puts $fp "plot '$datafile' $plot_instruction"
    puts $fp {yaxis_max = GPVAL_Y_MAX}
    puts $fp {ydata_max = GPVAL_DATA_Y_MAX}
    puts $fp {yaxis_min = GPVAL_Y_MIN}
    puts $fp {ydata_min = GPVAL_DATA_Y_MIN}
    puts $fp "if (yaxis_max < ($plot_padding * ydata_max)) \{"
    puts $fp "  set yrange \[(ydata_min / $plot_padding):($plot_padding * ydata_max)\] "
    puts $fp "  replot"
    puts $fp "\}"
    puts $fp "if (yaxis_min > (ydata_min / $plot_padding)) \{"
    puts $fp "  set yrange \[(ydata_min / $plot_padding):($plot_padding * ydata_max)\] "
    puts $fp "  replot"
    puts $fp "\}"
    puts $fp "set output '$scriptname.eps'"
    puts $fp "set terminal postscript eps color size 6in,4in"
    puts $fp "replot"
    puts $fp "set terminal wxt size 800,600"
    puts $fp "replot"
    close $fp
}

########################## Main entry point ##########################

set speaker_dict [dict create]

# Resonant frequency (Hz)
dict set speaker_dict fc 80.0

# DC resistance (ohms)
dict set speaker_dict re 7.6

# Mechanical Q
dict set speaker_dict qmc 8.0

# Electrical Q
dict set speaker_dict qec 0.73

# Starting frequency (Hz)
set fstart 20

# Ending frequency (Hz)
set fend 300

# Frequency points per decade
set perdecade 50

set freq_list [logrange $fstart $perdecade $fend]


set zcoil_mag_list [list]

foreach freq $freq_list {
    lappend zcoil_mag_list [magnum [zcoil $freq [dict get $speaker_dict re] \
					[dict get $speaker_dict fc] \
					[dict get $speaker_dict qmc] \
					[dict get $speaker_dict qec]]]
    
}


write_gnuplot_data data.dat $freq_list $zcoil_mag_list
write_gnuplot_script simimp.gp data.dat
# exec gnuplot -persist potplot.gp