/* -*- MaTX -*- * * NAME * nicholsp() - Plot the gain-phase response on the Nichols chart * * SYNOPSIS * nicholsp() * * nicholsp(G) * Rational G; * * nicholsp(G, w) * Rational G; * Array w; * * DESCRIPTION * nicholsp(G, w) plots the gain-phase response of SISO * continuous-time linear system on the Nichols chart. * * Vector w must contain the frequencies, in radians, at which the * Nichols response is to be evaluated. nicholsp() plots gain and phase * (in radian) with as many columns with length(w) rows. * * SEE ALSO * logspace, bode, and nyquist */ Func void nicholsp(G, w, ...) Rational G; Array w; { Integer i, j, Np, Nm, Na, win; Real phase_max, phase_min, gain_max, gain_min, ca, m, m2, al, th; Array G_gain, G_phase, G_Mp, G_Mm, G_M0; Array Mp, Mm, theta_p, theta_0, theta_m, Theta_p; Array Alpha, Theta_a, theta, Ga; List data_tits_p, data_cmds_p, data_tits_m, data_cmds_m; List data_tits_a, data_cmds_a; error(nargchk(1, 2, nargs, "nicholsp")); if (nargs < 2) { w = logspace(-2.0, 1.0); } if (1 <= nargs) { {G_gain, G_phase} = Bode_tf(G, w); G_gain = 20*log10(G_gain); gain_max = ceil(max(G_gain/5))*5; gain_min = floor(min(G_gain/5))*5; phase_max = (ceil(max(G_phase)/90)+1)*90; // deg phase_min = (floor(min(G_phase)/90)-1)*90; // deg if (gain_max < 40.0) { gain_max = 40.0; } if (-30.0 < gain_min) { gain_min = -30.0; } } else { gain_max = 40.0; gain_min = -30.0; } // Constant M loci // M = 0 theta_0 = linspace(-269.0, -91.0, 50); // deg G_M0 = Z(1,length(theta_0)); for (j = 1; j <= length(theta_0); j++) { ca = cos(theta_0(j)/180*PI); G_M0(1,j) = -1/(2*ca); } // M < 0 Mm = -[0.5, 1, 2, 4, 6, 10, 15, 20, 27, 33, 40]; // dB theta_m = linspace(-360.0, 0.0, 50); // deg Nm = length(Mm); G_Mm = Z(Nm,length(theta_m)); for (i = 1; i <= Nm; i++) { m = 10^(Mm(i)/20); m2 = m*m; for (j = 1; j <= length(theta_m); j++) { ca = cos(theta_m(j)/180*PI); G_Mm(i,j) = (m2/(m2 - 1))*(-ca - sqrt(ca*ca - 1 + 1/m2)); } } // M > 0 Mp = [12 6 4 3 2 1 0.5 0.25]; // dB theta_p = linspace(-270.0, -90.0); // deg theta_p = [theta_p fliplr(theta_p)]; // deg Np = length(Mp); G_Mp = Z(Np,length(theta_p)); Theta_p = Z(Np,length(theta_p)); for (i = 1; i <= Np; i++) { m = 10^(Mp(i)/20); m2 = m*m; for (j = 1; j <= length(theta_p); j++) { ca = cos(theta_p(j)/180*PI); if (j <= length(theta_p)/2) { if (ca*ca - 1 + 1/m2 >= 0.0) { Theta_p(i,j) = theta_p(j); G_Mp(i,j) = (m2/(m2 - 1))*(-ca + sqrt(ca*ca - 1 + 1/m2)); } else { if (j != 1) { Theta_p(i,j) = Theta_p(i,j-1); G_Mp(i,j) = G_Mp(i,j-1); } else { Theta_p(i,j) = 0.0; // temporary point G_Mp(i,j) = 0.0; // temporary point } } } else { if (ca*ca - 1 + 1/m2 >= 0.0) { Theta_p(i,j) = theta_p(j); G_Mp(i,j) = (m2/(m2 - 1))*(-ca - sqrt(ca*ca - 1 + 1/m2)); } else { Theta_p(i,j) = Theta_p(i,j-1); G_Mp(i,j) = G_Mp(i,j-1); } } } // Change the temporary points for (j = 1; j <= length(theta_p)/2; j++) { if (Theta_p(i,j) == 0.0 && G_Mp(i,j) == 0.0) { Theta_p(i,j) = Theta_p(i,length(theta_p)); G_Mp(i,j) = G_Mp(i,length(theta_p)); } } } // // Constant alpha loci // Alpha = [150, 120, 90, 60, 30, 20, 10, 5, 2]; Alpha = [-Alpha, fliplr(Alpha)]; theta = linspace(-360.0, 0.0); Na = length(Alpha); Ga = Z(Na,length(theta)); Theta_a = Z(Na, length(theta)); for (i = 1; i <= length(Alpha); i++) { al = Alpha(i)/180*PI; for (j = 1; j <= length(theta); j++) { th = theta(j)/180*PI; Ga(i,j) = (1/tan(al))*sin(th) - cos(th); Theta_a(i,j) = theta(j); if (Ga(i,j) < 0.0) { Ga(i,j) = EPS; Theta_a(i,j) = theta(j); } } } data_tits_p = makelist(Np); data_cmds_p = makelist(Np); for (i = 1; i <= Np; i++) { data_tits_p(i) = ""; data_cmds_p(i) = "2"; } data_tits_m = makelist(Nm); data_cmds_m = makelist(Nm); for (i = 1; i <= Nm; i++) { data_tits_m(i) = ""; data_cmds_m(i) = "2"; } data_tits_a = makelist(Na); data_cmds_a = makelist(Na); for (i = 1; i <= Na; i++) { data_tits_a(i) = ""; data_cmds_a(i) = "3"; } win = 3; mgplot_grid(win, 1); mgplot_xlabel(win, "phase [deg]"); mgplot_ylabel(win, "gain [dB]"); mgplot_cmd(win, "set xrange [-360:0]"); mgplot_cmd(win, "set xtics -360, 60, 0"); mgplot_cmd(win, "set yrange " + sprintf("[%g:%g]", gain_min, gain_max)); mgplot_text(win, "0dB", -100.0, 27.0, "center"); mgplot_text(win, "-0.5dB", -5.0, 26.0, "right"); mgplot_text(win, "-1dB", -5.0, 19.0, "right"); mgplot_text(win, "-2dB", -5.0, 13.0, "right"); mgplot_text(win, "-4dB", -5.0, 6.0, "right"); mgplot_text(win, "-6dB", -5.0, 1.5, "right"); mgplot_text(win, "-10dB", -5.0, -5.5, "right"); mgplot_text(win, "-15dB", -5.0, -12.0, "right"); mgplot_text(win, "-20dB", -5.0, -18.0, "right"); mgplot_text(win, "-27dB", -5.0, -25.0, "right"); mgplot_text(win, "0.25dB", -180.0, 32.0, "center"); mgplot_text(win, "0.5dB", -180.0, 26.0, "center"); mgplot_text(win, "1dB", -180.0, 20.0, "center"); mgplot_text(win, "2dB", -180.0, 15.0, "center"); mgplot_text(win, "3dB", -180.0, 11.5, "center"); mgplot_text(win, "4dB", -180.0, 9.5, "center"); mgplot_text(win, "6dB", -180.0, 7.0, "center"); mgplot_text(win, "12dB", -180.0, 3.5, "center"); mgplot_text(win, "-2deg", -70.0, 30.0, "center"); mgplot_text(win, "-5deg", -70.0, 22.0, "center"); mgplot_text(win, "-10deg", -70.0, 15.0, "center"); mgplot_text(win, "-20deg", -70.0, 9.0, "center"); mgplot_text(win, "-30deg", -70.0, 4.0, "center"); mgplot_text(win, "-60deg", -82.0, -3.0, "center"); mgplot_text(win, "-90deg", -100.0, -7.0, "center"); mgplot_text(win, "-120deg", -120.0,-10.0, "center"); mgplot_text(win, "-150deg", -145.0,-13.0, "center"); mgplot_text(win, "+2deg", -290.0, 30.0, "center"); mgplot_text(win, "+5deg", -290.0, 22.0, "center"); mgplot_text(win, "+10deg", -290.0, 15.0, "center"); mgplot_text(win, "+20deg", -290.0, 9.0, "center"); mgplot_text(win, "+30deg", -290.0, 4.0, "center"); mgplot_text(win, "+60deg", -278.0, -3.0, "center"); mgplot_text(win, "+90deg", -260.0, -7.0, "center"); mgplot_text(win, "+120deg", -240.0,-10.0, "center"); mgplot_text(win, "+150deg", -215.0,-13.0, "center"); if (1 <= nargs) { mgplot(win, G_phase, G_gain, {"G(s)"}, {"1"}); } mgreplot(win, [0], [0], {"alpha loci"}, {"3"}); mgreplot(win, theta_0, 20*log10(abs(G_M0)), {"M loci"}, {"2"}); mgreplot(win, theta_m, 20*log10(abs(G_Mm)), data_tits_m, data_cmds_m); mgreplot(win, Theta_p, 20*log10(abs(G_Mp)), data_tits_p, data_cmds_p); mgreplot(win, Theta_a, 20*log10(abs(Ga)), data_tits_a, data_cmds_a); }