steps_per_dec = 6;
decades = 6;
start_freq = 0.1;

// The transfer function
function foo=G(w)
	D = %i * w;
	foo = (D + 5) / (D^2 + 100*D + 10000);
endfunction

fd = mopen("data.txt", "w");
for step = 0:(steps_per_dec * decades),
	f = start_freq * 10 ^ (step/steps_per_dec);	// calculate the next frequency
	w = f * (2 * %pi);		// convert the frequency to radians
	[gain, phase] = polar(G(w));	// get the result and convert to magnitude/angle
	gaindb = 20*log10(gain);	// convert to dB
	phasedeg = 180 * phase / %pi;	// convert to degrees
	mfprintf(fd, "%f, %f, %f \n", f, gaindb, phasedeg);
end
mclose(fd);


// To graph it directly
D=poly(0,'D');
h=syslin('c', (D + 5) / (D^2 + 100*D + 10000) );
bode(h, 0.1, 1000, 'Sample Transfer Function');