# Initial boilerplate for setting up gnuplot files for impedance calculations

j=sqrt(-1)

zc(f,c) = 1/(j*2*pi*f*c)	# capacitors
zl(f,l) = j*2*pi*f*l		# inductors

# Weird, non-linear "semi-impedance" for modeling loudspeakers
zsemi(f,l,alpha) = j * (2*pi*f)**alpha *l

zpar(z1,z2) = z1*z2/(z1+z2)  # parallel impedances

divider(zup,zdown) = zdown/(zup+zdown) # voltage divider gain
noninvert(zfeedback,zgnd) = 1 + zfeedback/zgnd	# noninverting amplifier gain
invert(zfeedback,zin) = - zfeedback/zin	# inverting amplifier gain

set logscale xy
set title "Impedance---CHANGE THIS TITLE"
set ylabel "magnitude of impedance (ohms)"
set xlabel "frequency (Hz)"
set samples 1000	# plot a smooth curve for the model
			# even if it has sharp spikes
set key bottom left	# find a place to stick the key

# SELECT A FREQUENCY RANGE OF INTEREST
set xrange [1:50e3]	# for audio signals
set xrange [1:2e6]	# for full range of multimeter (plus a bit)
set xrange [0.1:100]	# for pulse measurements

# MAKE FUNCTIONS FOR YOUR MODELS, WITH APPROPRIATE PARAMETERS
zmodel(f,r,c) = zpar(r, zc(f,c))	# r || c

# REMEMBER THAT WE ARE MAINLY MEASURING OR FITTING MAGNITUDE OF IMPEDANCE.
# YOU CAN MAKE SOME OF THE PARAMETERS BE GLOBAL VARIABLES:

azmodel(f) = abs(zmodel(f,R1,C1))

# SET VALUES FOR YOUR PARAMETERS (EITHER COMPONENT VALUES OR INITIAL GUESSES)

R1= 5.6e6	# 5.6Mohm
C1= 1e-9	# 0.001 uF

# FIT YOUR MODEL TO DATA (IF ANY), MAKING SURE THAT THE 
# CALCULATION FOR THE DATA RESULTS IN AN IMPEDANCE,

# USE LINEAR FIT TO EMPHASIZE HIGH IMPEDANCE REGION
# fit azmodel(x) "data_file" using 1:($3/$2*47) via R1,C1

# USE LOG FIT TO EMPHASIZE LOW IMPEDANCE REGION
# fit log(azmodel(x)) "data_file" using 1:(log($3/$2*47)) via R1,C1

# MAKE A NAME FOR THE MODEL (ECHOING PARAMETERS)
model_name = sprintf("%.3e ohm || %.3e F",R1,C1)


unset arrow	# turn off any stray arrows or lines
unset label	# turn off any stray labels
set xtics	# use default xtics

# ADD LABELS, ARROWS, OR XTICS AS NEEDED

corner = 1./(2*pi*R1*C1)
set arrow from corner,1e6 to corner, 1e7 nohead
set label sprintf("corner=%.2f Hz",corner) at corner,3e6 right

z60 = azmodel(60) # magnitude of impedance at 60Hz
set label sprintf("%.3e ohms at 60Hz", z60) at 60,z60 right

plot azmodel(x) title model_name