/* check various formulas for a damped harmonic oscillator including area under the lineshape (based on the power, i.e. gain *squared*) as well as peak position, FWHM, and a few other things */ /* invoke with either: batch("rlc.max"); maxima -b rlc.max *********/ Gsq : 1 / (Q^-2 * phi^2 + (phi^2 - 1)^2); assume (Q > 1/2); /* required for integral */ area_of_Gsq : defint(Gsq, phi, 0, inf); G2peak : 4 * Q^4 / (4*Q^2 - 1); Gpeak : sqrt(G2peak); /* location of peak */ Psipeak : (2*Q^2 - 1) / (2 * Q^2); if (ratsimp(ev(Gsq, phi=sqrt(Psipeak)) - G2peak) # 0) then break("wrong peak gain")$ /* now, for the bandwidth: */ Phi1 : ratsimp(sqrt(Psipeak - 1/Gpeak)); Phi2 : ratsimp(sqrt(Psipeak + 1/Gpeak)); if (ratsimp(ev(Gsq/G2peak, phi=Phi1)) # 1/2) then break("not lower halfway")$ if (ratsimp(ev(Gsq/G2peak, phi=Phi2)) # 1/2) then break("not upper halfway")$ /* complicated and ugly */ FWHM : ratsimp(Phi2 - Phi1); /* convert to 1/Q for Taylor series: */ FWHMbar : ev(FWHM, Q=1/Qbar)$ taylor(FWHMbar, Qbar, 0, 5); /* estimate range of validity */ /* by finding where the lower root goes to zero */ /* product of the roots is zero iff the lower root is zero */ zipx : ratsimp(Psipeak^2 - 1/G2peak)$ zip : zipx * 2 * Q^4; /* the smallest possible Q that has any real FWHM */ /* found by finding the roots of zip */ Qmin : sqrt(1 + sqrt(1/2)); /* verify that it is a root */ if (ratsimp(ev(zip, Q=Qmin)) # 0) then break("not a root")$ float(Qmin); float(ev(FWHM, Q=10)); float(ev(FWHM, Q=2)); float(ev(FWHM, Q=Qmin)); float(ratsimp(ev(FWHM, Q=13/10))); /* slightly too small */ /* location of the actual peak of the resonance */ phi_res : sqrt(1 - 1/2/Q^2)$ /* another sanity check */ /* at Qmin, the peak height should be 2 */ /* that's why the concept of *half* max becomes invalid */ simp1 : ratsimp(ev(Gsq, phi=phi_res))$ if (ratsimp(ev(simp1, Q=Qmin)) # 2) then break("wrong low peak")$ /* Area under the curve of plain voltage gain (not gain squared). Doesn't work. Maxima doesn't know how to integrate it. Not important. All we really care about is gain squared anyway. xx gplain : sqrt(Gsq); xx area_plain : defint(gplain, phi, 0, inf); */ peak : ratsimp(ev(Gsq, phi = phi_res)); float(ev(phi_res, Q=1)); float(ev(peak, Q=1)); float(ev(Gsq,phi=.707,Q=1)); /* resden is resonant denominator, almost in Lorentz form */ resden : psi^2 - (2-1/Q^2)*psi + 1; /* inverse of G2peak, from completing the square in resden */ Gm2pb : resden - (psi - Psipeak)^2; /* check that Lorentz form is correct */ if (ratsimp(Gm2pb * G2peak) # 1) then break("bad inverse G2")$ /* check alternate form for Psipeak */ if (ratsimp((1-1/2/Q^2) - Psipeak) # 0) then break("Psipeak not the same")$ /* All ok if we make it to here. */ "OK"$ quit()$