Q = 0.25;
L = 5;
[D1, fn] = ndgrid(0.1:0.1:0.5, 0.5:0.1:1.5);
D1_fn = D1./fn; %optimization
a = 1 - cos(2*pi*D1_fn);
x = (pi*D1) ./ (fn*L) .* sin(2*pi*D1_fn);
y = cos(2*pi*D1_fn) .* ((4*Q/pi*fn) + ((pi^2*D1_fn./fn*L) .* (0.5-D1)) - 1);
z = (pi^2*D1_fn./fn*L) .* (0.5-D1);
m = 4*Q/pi * fn;
deno = 1+x+y+z+m;
Gdc1 = a ./ deno;
surf(D1, fn, Gdc1, 'edgecolor', 'none')
