Snippets
Created by
John Peck
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 | #!/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
|
Comments (0)
You can clone a snippet to your computer for local editing. Learn more.