Where did the LP to BP transformation $ \frac{\omega ^{2}-\omega_0^{2}}{\omega\cdot BW} $ come from?

Question

I'm not sure this question is appropriate for this site, but I feel that if anyone knows about a filter, a radio ham knows! And as a ham myself I feel I need to understand this important area better. To that end, I've been reading various filter texts (Zverev and Dimopoulos mainly).

I'm at the point where the low-pass approximation has been developed (I'm going to chose a Butterworth for this question) and I'm moving to the transformation of a low-pass to other types (high-pass, band-pass, etc).

Lets take a Butterworth low-pass given in the book by Dimopoulos: $$ \frac{17.5514}{s^3+3.2743\cdot s^2+5.3604 \cdot s+4.3878}$$

If you have Maxima installed, you can try it yourself:

load(draw);

maxima_tempdir:"C:/maxima-5.41.0a";

H_lpbut(s):=17.5514/(s^3+3.2743*s^2+5.3604*s+4.3878);

G_lp(omega):=abs(subst([s=%i*omega],H_lpbut(s)));

Omega_3dB(N):=(1/beta_m())^(1/N);

draw2d(
        xlabel="Omega" 
        ,ylabel="Lowpass G(Omega)"
        ,xaxis=true
        ,yaxis=true
        ,xrange=[-10,10]
        ,yrange=[0,5]
        ,explicit(G_lp(omega),
                    omega,-10,10
                    )
        ,color=red
        ,points([[Omega_3dB(3), G_lp(Omega_3dB(3))]])
        ,points_joined=true
        ,line_type=dashes
        ,points([
                [0,G_lp(Omega_3dB(3))]
                ,[Omega_3dB(3), G_lp(Omega_3dB(3))]
                ,[Omega_3dB(3), 0]
            ])
    );

As you can see, I get the magnitude of the transfer function and plot from the normalized frequencies -10 to 10. Notice how wonderfully symmetrical the plot is.

That sort of makes sense, we're plugging in negative frequencies (from -10 to zero) and we end up with a mirror image of the low-pass giving a band-pass shape.

In fact, that's pretty much how the lp to bp transform works, only this time they use a different function:

$$ \frac{\omega ^{2}-\omega_0^{2}}{\omega\cdot BW} $$

Again, using Maxima we can explore how this function maps band-pass frequencies to normalized low-pass frequencies. I've taken some liberties here to highlight it's properties (it uses a log scale and negative infinity is now positive infinity - but you get the idea). The function goes to infinity at one end (omega=0) and takes a more leisurely trip to infinity at the other (omega=inf). It converges at zero (or bp centre frequency), which is exactly what required for bp transformation - so I get that bit.

bpt(omega,omega_0,BW):=(omega^2-omega_0^2)/(omega*BW);

draw2d(
        xrange=[0,1000]
        ,yrange=[0.0001,1000]
        ,logy=true
        ,explicit(sqrt(bpt(2*%pi*f,2*%pi*500,2*%pi*100)^2),
                    f,0,1000
                    )
    );

and finally, here's the transformed band-pass using the standard lp-to-bp transform:

Notice it isn't as symmetrical as the previous plot (using -10 to 10).

But, here's what I don't get...why would they use a function like this, which is obviously non-symmetrical, to produce the band-pass? Why not use a function that works to produce the nice linear -10 to 10 input that we saw earlier?