¿Cómo calculo la impedancia de un condensador (bypass) (para simulación)?

Quería construir una simulación de condensadores de derivación en Matlab, especialmente para jugar con pares de condensadores de diferentes valores. Así que pensé en tomar un modelo de un capacitor real como una conexión en serie de capacitor con una resistencia (ESR) y un inductor (ESL) y derivar la fórmula

$$Z=\mathrm{ESR}+j2\pi f\cdot\mathrm{ESL}+\frac{1}{j2\pi fC}$$ separando en parte real e imaginaria: $$Z=\mathrm{ESR}+j\left(2\pi f\cdot\mathrm{ESL}-\frac{1}{2\pi fC}\right)$$ ahora toma la magnitud $$|Z|=\sqrt{\mathrm{ESR}^2+\left(2\pi f\cdot\mathrm{ESL}-\frac{1}{2\pi fC}\right)^2}$$

Ahora programé el Matlab. primero la función que calcula la impedancia de la tapa a la frecuencia o frecuencias dadas

function [ Z ] = calcCapImp( f,R, L,C )  
Z=sqrt(R^2+(2*pi*f*L-1./(2*pi*C*f)).^2);
end

El código para calcular y graficar la impedancia de una combinación de dos capacitores en Ohms y dBs:

clear all;
f=10000:1000:1e9;
Z=(zeros(size(f)));
for i=1:size(f,2)
Z1=calcCapImp(f(i),10e-3,700e-12,100e-9);%calculate first cap
Z2=calcCapImp(f(i),10e-3,1250e-12,4.7e-6);%calculate second  cap
Z(i)=(Z1*Z2)/(Z1+Z2);%Total impedance is parallel combination
end
loglog(f,Z);
ylabel('Resistance, Ohm');
xlabel('Frequency, Hz');
grid on;
figure;
Z=20*log10(Z);%convert to decibels
semilogx(f,Z);
xlabel('Frequency, Hz');
ylabel('Resistance, dB');
grid on;

Si bien los resultados parecen plausibles,

No puedo reproducir el gráfico que he visto en otros lugares.

En particular, el código anterior intenta recrear la Figura 2 desde aquí (mismo rango):

No obtengo un pico tan dramático a ~ 10 MHz como en la figura externa mencionada anteriormente (aproximadamente 10 veces/20 dB más pequeña)

¿He hecho bien los cálculos? ¿Hay algún error en mi simulación?

EDITAR: aquí el código y el gráfico corregidos según la respuesta de The Photon a continuación

function [ Z ] = calcCapImp( f,R, L,C )

Z=R+1i*(2*pi*f*L-1./(2*pi*C*f));

end

clear all;
f=10000:1000:1e9;
Z=(zeros(size(f)));
for i=1:size(f,2)
Z1=calcCapImp(f(i),10e-3,700e-12,100e-9);%calculate first cap
Z2=calcCapImp(f(i),10e-3,1250e-12,4.7e-6);%calculate second  cap
Z(i)=abs((Z1*Z2)/(Z1+Z2));%Total impedance is parallel combination, take the magnitude
end

loglog(f,Z);
ylabel('Impedance, Ohm');
xlabel('Frequency, Hz');
title('Impedance in Ohms');
grid on;
figure;
Zdb=20*log10(Z);%convert to decibels
semilogx(f,Zdb);
xlabel('Frequency, Hz');
ylabel('Impedance, dB');
title('Impedance in decibels');
grid on;

Ahora el resultado es prácticamente idéntico a los utilizados como referencia.