set logscale xy
set xlabel "frequency [Hz]"
set xrange[10:1e6]
set key bottom left Left
set samples 1000
set grid

j=sqrt(-1)
conjugate(a) = real(a) - j*imag(a)
phase(v) = imag(log(v))
real_power(V,Z)= real(V*conjugate(V/Z)) # V is RMS voltage phasor

divider(Zup,Zdown) = Zdown/(Zup+Zdown)
Zc(f,C) = 1/(j*2*pi*f*C)  # f= frequency [Hz], C=cap [farads]
Zl(f,L) = j*2*pi*f*L
Zpar(z1,z2) = z1*z2/(z1+z2)

#model for impedance of loudspeaker
Rnom=7.65
Lhi= 5.78e-3    ;    Lhipow = -0.376
Ls=2.91e-3      ;    Cs=389e-6   ;  Rs=15.5

#impedance of loudspeaker
Zloud(f) = Rnom+Zl(f,Lhi*f**Lhipow)+Zpar(Rs,Zpar(Zl(f,Ls),Zc(f,Cs)))
Rseries = 0.0 # extra resistor for displaying current on scope

# voltage divider with Lfilt and Cfilt
gain(f,Lfilt,Rfilt,Cfilt) = divider(Zl(f,Lfilt)+Rfilt, Zpar(Zloud(f)+Rseries,Zc(f,Cfilt)))

# total impedance, in case we want to plot that also
Ztotal(f,Lfilt,Rfilt,Cfilt) = Zl(f,Lfilt)+Rfilt + Zpar(Zloud(f)+Rseries,Zc(f,Cfilt))

# real power into loudspeaker with voltage from amplifier V at frequency f
loud_power(V,f,Lfilt,Rfilt,Cfilt) = real_power(V*gain(f,Lfilt,Rfilt,Cfilt)*divider(Rseries,Zloud(f)), Zloud(f))
 
nofilter_power(V,f)=real_power(V*divider(Rseries,Zloud(f)), Zloud(f))
 
V_supply=3.3    # power supply
V_RMS = V_supply/sqrt(2)        # RMS voltage for largest sine wave from supply
Lfilt=220e-6;   Rfilt=0.252     # parameters for AIUR-06-221 inductor
# Lfilt=270e-6;   Rfilt=0.270     # parameters for AIUR-06-271 inductor
# Lfilt=330e-6;   Rfilt=0.394     # parameters for AIUR-06-331 inductor

set title sprintf("Low-pass LC filter for loudspeaker with %.3gohm series resistor",Rseries)

half_power=0.5*loud_power(V_RMS,300,Lfilt,Rfilt,3.3e-6) # half-power point

set yrange [1e-4:40]
set ylabel sprintf('Real power to speaker with +- %.3gV supply [W]',V_supply)
plot half_power title sprintf("half-power=%.3gW",half_power), \
        nofilter_power(V_RMS,x) title "no filter", \
        loud_power(V_RMS,x,Lfilt,Rfilt,0.1e-6) title sprintf("L=%.3guH C=%.3guF", Lfilt*1e6, 0.1), \
        loud_power(V_RMS,x,Lfilt,Rfilt,1e-6) title sprintf("L=%.3guH C=%.3guF", Lfilt*1e6, 1), \
        loud_power(V_RMS,x,Lfilt,Rfilt,2.2e-6) title sprintf("L=%.3guH C=%.3guF", Lfilt*1e6, 2.2), \
        loud_power(V_RMS,x,Lfilt,Rfilt,3.3e-6) title sprintf("L=%.3guH C=%.3guF", Lfilt*1e6, 3.3), \
        loud_power(V_RMS,x,Lfilt,Rfilt,4.7e-6) title sprintf("L=%.3guH C=%.3guF", Lfilt*1e6, 4.7), \
        loud_power(V_RMS,x,Lfilt,Rfilt,10e-6) title sprintf("L=%.3guH C=%.3guF", Lfilt*1e6, 10)