Listing of GNU/Octave code for RLC resonator. Used to generated the graphs in my RLC resonator article. Supplements the article RLC resonator.
Appendix
Bode magnitude, in GNU/Octave
GNU/Octave code for RLC resonator clc; close all; clear all; format short eng L=47e-3; # 47mH C=47e-9; # 47nF #Rvector = [3.9e3]; # separate real poles #Rvector = [2e3]; # coinciding real poles Rvector = [18 220 820]; # conjugate complex poles f=logspace(1,6,200); w=2*pi*f; for R = Rvector wn=1/sqrt(L*C); zeta=(R/2)*sqrt(C/L) if (zeta<1 ) # complex conjugate poles on left side of s-plane u=20*log10(R/L) + 20*log10(w) - 20*log10(sqrt((wn.^2-w.^2).^2+(2*zeta*wn*w).^2)); hold on; h=semilogx(f,u); hold off; hold on; hold off poles = [-wn*zeta+sqrt(1-zeta^2)*j, -wn*zeta-sqrt(1-zeta^2)*j ]; endif if (zeta == 1) # coinciding real poles p=-R/(2*L); #1/sqrt(L*C); u=20*log10(R/L) + 20*log10(w) - 40*log10(sqrt(w.^2+p.^2)); hold on; h=semilogx(f,u); plot([wn/(2*pi) wn/(2*pi)], get(gca,'YLim'),'k--'); text(wn/(2*pi),5,'|p|/2\pi'); f1=-p/(2*pi); fmin=min(f); fmax=max(f); #asymp1=0 - 20*log10((f1-fmin)/fmin); #asymp2=0 - 20*log10((fmax-f1)/f1); #plot([fmin f1 fmax],[asymp1 0 asymp2],'k--'); hold off poles=[-wn*zeta -wn*zeta]; endif if (zeta>1) # separate real poles p1=wn*(-zeta+sqrt(zeta^2-1)); p2=wn*(-zeta-sqrt(zeta^2-1)); u=20*log10(R/L) + 20*log10(w) – 20*log10(sqrt(w.^2+p1.^2)) – 20*log10(sqrt(w.^2+p2.^2)); figure(1); hold on; h=semilogx(f,u); hold off; hold on; f1=-p1/(2*pi); f2=-p2/(2*pi); plot([f1 f1], get(gca,’YLim’),’k–‘); plot([f2 f2], get(gca,’YLim’),’k–‘); fmin=min(f); fmax=max(f); asymp1=0 – 20*log10((f1-fmin)/fmin); asymp2=0 – 20*log10((fmax-f2)/f2); plot([min(f) f1 f2 fmax],[asymp1 0 0 asymp2],’k–‘); text(f1,5,’|p1|/2\pi’); text(f2,5,’|p2|/2\pi’); hold off poles=[p1+0j p2+0j]; endif #figure(2); #grid off; ##axis([-1e4 1e4 -1e4 1e4]); ##hold on; #zplane([],poles’); ##hold off; endfor figure(1); grid off; axis([min(f) max(f) -100 20]); xlabel(‘frequency [Hz]’); ylabel(’20log| H(t)|’); leg=strread(num2str(Rvector,4),’%s’); if (zeta>=1) leg=[leg;’asymptote’]; endif t=[‘Bode Magnitude in dB(f’]; t2=[‘), C=’ num2str(C*1e9) ‘\muF, L=’, num2str(L*1e3),’mH’]; if(length(Rvector)==1) t=[t t2 ‘, R=’ num2str(R/1e3) ‘k\Omega’] else t=[t ‘,R’ t2]; legend(leg); endif if (zeta<1) hold on; plot([wn/(2*pi) wn/(2*pi)], get(gca,'YLim'),'k--'); text(wn/(2*pi),5,'\omega_n/2\pi'); hold off; endif t = [t ]; title(t, "fontsize";, 15);[/code]