#!/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: $scriptname \[options] filename"
set options {
    {t "Plot total impedance"}
    {r "Plot Rec"}
    {p "Plot total impedance phase"}
}

# 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} {message optdict} {
    # message will be the human-readable error.  optdict will be the
    # options dictionary present when the command was attempted.
    puts $message
    exit 1
}

# After cmdline is done, argv will point to the last argument
if {[llength $argv] == 1} {
    set infile $argv
} else {
    # Input file hasn't been specified.
    puts "No input file specified"
    puts $usage
    exit
}



# Constants
array set constants {
    pi 3.1415926535897932385
}

########################## 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


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

proc deg_to_rad {degrees} {
    # Return the angle in radians
    global constants
    set radians [expr $degrees * ($constants(pi) / 180)]
    return $radians
}

proc rad_to_deg {radians} {
    # Return the angle in degrees
    global constants
    set degrees [expr $radians * (180 / $constants(pi))]
    return $degrees
}

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

proc get_raw_dict {fp} {
    # Return a dictionary of dictionaries.  Each key in the topmost
    # dictionary will be a frequency point.
    #
    # Arguments:
    #
    # fp -- (open) I/O channel to raw CGR datafile
    set line_list [split [read $fp] "\n"]
    set raw_dict [dict create]
    # Skip first heading line
    foreach line [lrange $line_list 1 end] {
	set param_list [split $line ","]
	if {[llength $param_list] == 0} {
	    # We've reached the end
	    break
	}
	set frequency [lindex $param_list 0]
	set xfer_db [lindex $param_list 1]
	set phase_deg [lindex $param_list 2]
	set quality [lindex $param_list 3]
	# Create a frequency key with a parameter dictionary
	dict set raw_dict $frequency [dict create]
	# Get a handle to this new dictionary
	set freq_dict [dict get $raw_dict $frequency]
	set xfer_rad [deg_to_rad $phase_deg]
	set xfer_mag [db_to_ratio $xfer_db]
	set xfer_cpx [complex [expr $xfer_mag * cos($xfer_rad)] \
			  [expr $xfer_mag * sin($xfer_rad)]]
	dict set freq_dict "xfer_cpx" $xfer_cpx
	dict set freq_dict "phase_deg" $phase_deg
	dict set freq_dict "quality" $quality
	# Save the dictionary back to the measurement dictionary
	dict set raw_dict $frequency $freq_dict              
    }
    return $raw_dict
}

proc get_z_mag_list {raw_dict} {
    # Return a list of impedance magnitudes in ohms
    #
    # Arguments:
    #
    # raw_dict -- Dictionary created by get_raw_dict
    global cal_dict
    set zmag_list [list]
    foreach frequency [dict keys $raw_dict] {
	set freq_dict [dict get $raw_dict $frequency]
	set xfer_cpx [dict get $freq_dict xfer_cpx]
	set xfer_inv_cpx [/ [complex 1 0] $xfer_cpx]
	set rsense_cpx [complex [dict get $cal_dict rsense] 0]
	set z_norm_cpx [- $xfer_inv_cpx [complex 1 0]]
	set z_cpx [* $rsense_cpx $z_norm_cpx]
	set z_mag [magnum $z_cpx]
	lappend z_mag_list $z_mag
    }
    return $z_mag_list
}

proc get_r_mag_list {raw_dict} {
    # Return a list of resistance magnitudes in ohms
    #
    # Arguments:
    #
    # raw_dict -- Dictionary created by get_raw_dict
    global cal_dict
    global seed_dict
    set rmag_list [list]
    foreach frequency [dict keys $raw_dict] {
	set freq_dict [dict get $raw_dict $frequency]
	set xfer_cpx [dict get $freq_dict xfer_cpx]
	set xfer_inv_cpx [/ [complex 1 0] $xfer_cpx]
	set rsense_cpx [complex [dict get $cal_dict rsense] 0]
	set re_cpx [complex [dict get $seed_dict re] 0]
	set z_norm_cpx [- $xfer_inv_cpx [complex 1 0]]
	set z_cpx [* $rsense_cpx $z_norm_cpx]
	# Subtract off Re to leave the parallel rlc combination
	set rlc_cpx [- $z_cpx $re_cpx]
	set y_cpx [/ [complex 1 0] $rlc_cpx]
	set r_inv [real $y_cpx]
	set r_mag [expr double(1)/$r_inv]
	lappend r_mag_list $r_mag
    }
    return $r_mag_list
}

proc get_phase_list {raw_dict} {
    # Return a list of phase angles in degress
    #
    # Arguments:
    #
    # raw_dict -- Dictionary created by get_raw_dict
    global cal_dict
    set phase_list [list]
    foreach frequency [dict keys $raw_dict] {	
	set freq_dict [dict get $raw_dict $frequency]
	set xfer_cpx [dict get $freq_dict xfer_cpx]
	set xfer_inv_cpx [/ [complex 1 0] $xfer_cpx]
	set rsense_cpx [complex [dict get $cal_dict rsense] 0]
	set z_norm_cpx [- $xfer_inv_cpx [complex 1 0]]
	set z_cpx [* $rsense_cpx $z_norm_cpx]
	set phase_rad [arg $z_cpx]
	lappend phase_list [rad_to_deg $phase_rad]
    }
    return $phase_list
}

proc db_to_ratio {dbval} {
    # Return the fractional ratio value corresponding to the decibel
    # power value
    #
    # Arguments:
    #
    # dbval -- 20 * log10 (ratio)
    set ratio [expr 10 ** (double($dbval)/20)]
    return $ratio
}

proc write_gnuplot_data {filename frequency_list impedance_list} {
    # Write a data file that can be plotted by gnuplot
    #
    # Arguments:
    #
    # filename -- Name of the output data file
    # frequency_list -- List of measurement frequencies
    # impedance_list -- List of impedances corresponding to frequencies
    global params
    # Open the output file
    try {
	set fp [open $filename w]
    } trap {} {message optdict} {
	puts $message
	exit
    }
    foreach frequency $frequency_list impedance $impedance_list {
	puts $fp "$frequency $impedance"
    }
    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 seed_dict
    global params
    # Open the output file
    try {
	set fp [open $filename w]
    } trap {} {message optdict} {
	puts $message
	exit
    }
    # 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
    puts $fp "reset"
    puts $fp "set terminal $gnuplot_terminal size 800,600"
    puts $fp "set format y '%0.0s %c'"
    puts $fp "set logscale x"
    puts $fp "set format x '%0.0s %c'"
    puts $fp "set xlabel 'Frequency (Hz)'"
    if $params(r) {
	puts $fp "set ylabel 'Rec (ohms)'"
    } elseif $params(p) {
	puts $fp "set ylabel 'Phase (deg)'"
    } else {
	puts $fp "set ylabel '|Z| (ohms)'"	
    }
    puts $fp "set nokey"
    puts $fp "set title '$scriptname'"
    if $params(r) {
	
    } elseif $params(p) {
	
    } else {
	# The voice coil impedance curve for fitting
	puts $fp "re = [dict get $seed_dict re]"
	puts $fp "fc = [dict get $seed_dict fc]"
	puts $fp "qec = [dict get $seed_dict qec]"
	set outstr "fitcurve(x, re, fc, qmc, qec) = re * abs( "
	append outstr "(-x**2 + {0,1}*fc*(1/qmc + 1/qec)*x + fc**2) / "
	append outstr "(-x**2 + {0,1}*(fc/qmc)*x + fc**2))"
	puts $fp $outstr
	set outstr "fit fitcurve(x, re, fc, qmc, qec) '$datafile' using 1:2 "
	append outstr "via fc, qmc, qec"
	puts $fp $outstr	
    }
    set plot_instruction "using 1:2 with lines"
    puts $fp "plot '$datafile' $plot_instruction"
    if $params(r) {
	
    } elseif $params(p) {
	
    } else {
	puts $fp "replot fitcurve(x, re, fc, qmc, qec) with lines"	
    }
    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 "\}"
    # Start annotations
    set anno_y 0.5
    set anno_x 0.7
    set anno_incr 0.05
    if $params(r) {
	
    } elseif $params(p) {
	
    } else {
	puts $fp "set label sprintf('Re = %0.1f Ohms',re) at graph $anno_x,$anno_y"
	set anno_y [expr $anno_y - $anno_incr]
	puts $fp "set label sprintf('fc = %0.0f Hz',fc) at graph $anno_x,$anno_y"
	set anno_y [expr $anno_y - $anno_incr]
	puts $fp "set label sprintf('Qec = %0.2f',qec) at graph $anno_x,$anno_y"
	set anno_y [expr $anno_y - $anno_incr]
	puts $fp "set label sprintf('Qmc = %0.2f',qmc) at graph $anno_x,$anno_y"	
    }

    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
}

set cal_dict [dict create]
# Reference resistor in ohms
dict set cal_dict rsense 10

set seed_dict [dict create]

# DC resistance of driver.  You should actually know this.
dict set seed_dict re 7.8

# Resonant frequency
dict set seed_dict fc 100

# Electrical Q
dict set seed_dict qec 0.8

# Mechanical Q
dict set seed_dict qmc 8.0

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

# Open the input file
try {
    set fp [open $infile r]
} trap {} {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
}

# Get the measurement dictionary
set raw_dict [get_raw_dict $fp]
set freq_list [dict keys $raw_dict]
if $params(r) {
    # We're calculating resistance
    set r_mag_list [get_r_mag_list $raw_dict]
    write_gnuplot_data "resistance.dat" $freq_list $r_mag_list
    write_gnuplot_script "plimp.gp" "resistance.dat"
    
} elseif $params(p) {
    # We're plotting phase
    set phase_list [get_phase_list $raw_dict]
    write_gnuplot_data "phase.dat" $freq_list $phase_list
    write_gnuplot_script "plimp.gp" "phase.dat"
} else {
    # We're looking at total impedance
    set z_mag_list [get_z_mag_list $raw_dict]
    write_gnuplot_data "linmag.dat" $freq_list $z_mag_list
    write_gnuplot_script "plimp.gp" "linmag.dat"
}


close $fp