% bodepowers.m - Bode plots of powers of k
% RMM, 24 Oct 05 (based on bodederint.m)

aminit;					% initialize MATALB

% Set up a function of the right form
g='s.^k';

% Set up the space of frequencies to use
w=logspace(-1,1,500); s=i*w;

% Evaluate the different powers of s
k=-2; gv=eval(g); gainm2=abs(gv); phasem2=180/pi*unwrap(angle(gv)) - 360;
k=-1; gv=eval(g); gainm1=abs(gv); phasem1=180/pi*unwrap(angle(gv));
k=0;  gv=eval(g); gainze=abs(gv); phaseze=180/pi*unwrap(angle(gv));
k=1;  gv=eval(g); gainp1=abs(gv); phasep1=180/pi*unwrap(angle(gv));
k=2;  gv=eval(g); gainp2=abs(gv); phasep2=180/pi*unwrap(angle(gv));

% Bode plots for negative powers of s
figure(1); clf; axes('position', [0.13 0.4096 0.3347 0.5]);
[magh, phaseh] = ambode(tf([1], [1 0 0]));	% Start with a double integrator

% Plot the gain
subplot(magh); hold on;
loglog(w, gainze, 'Linewidth', AM_bode_linewidth, 'Color', 'b');
loglog(w, gainm1, 'Linewidth', AM_bode_linewidth, 'Color', [0 0.5 0]);
loglog(w, gainm2, 'Linewidth', AM_bode_linewidth, 'Color', 'r');
amaxis([0.1 10 0.01 100]);
amtext(0.25, 40, 'sm2');
amtext(0.12, 16, 'sm1');
amtext(0.12, 0.5, 'sze');

% Plot the phase (note change of sign on phasep2, to get -180)
subplot(phaseh); hold on;
semilogx(w, phaseze, w, phasem1, w, phasem2, 'Linewidth', AM_bode_linewidth);
set(gca,'Ytick',[-180 0 180]);
amaxis([0.1 10 -200 200]);
amtext(0.15, -180 + 35, 'sm2');
amtext(0.15, -90 + 35, 'sm1');
amtext(0.15, 0 + 35, 'sze');

amprint('bodepowers-neg.eps');

% Bode plots for positive powers of s
figure(2); clf; axes('position', [0.13 0.4096 0.3347 0.5]);
[magh, phaseh] = ambode(tf([1 0 0], [1]));	% Start with a double integrator

subplot(magh); hold on;
loglog(w, gainze, 'Linewidth', AM_bode_linewidth, 'Color', 'b');
loglog(w, gainp1, 'Linewidth', AM_bode_linewidth, 'Color', [0 0.5 0]);
loglog(w, gainp2, 'Linewidth', AM_bode_linewidth, 'Color', 'r');
amaxis([0.1 10 0.01 100]);
amtext(0.12, 2, 'sze');
amtext(0.12, 0.3, 'sp1');
amtext(0.2, 0.03, 'sp2');

subplot(phaseh); hold on;
pl=semilogx(w, phaseze, w, phasep1, w, phasep2, 'Linewidth', AM_bode_linewidth);
set(gca,'Ytick',[-180 0 180]);
amaxis([0.1 10 -200 200]);
amtext(0.15, 90 -25, 'sp1');
amtext(0.15, 180 - 25, 'sp2');
amtext(0.15, 0 - 25, 'sze');

amprint('bodepowers-pos.eps');
