This is chap6.m in view mode; [Download] [Up]
(*
**********************************************************
COSY_PAK package chap6.m
**********************************************************
*)
If [ TrueQ[ $VersionNumber >= 2.0 ],
Off[ General::spell ];
Off[ General::spell1 ] ]
(* B E G I N P A C K A G E *)
BeginPackage["COSYPAK`chap6`","Algebra`ReIm`",
"Graphics`Graphics`"]
MagPlot::usage =
"MagPlot[Transf, s, {w,wmin,wmax}, gopts]: Plots the magnitude part of Bode plot of transfer function Transf(s) in decibel(dB), from frequency w=wmin to w=wmax, both > 0 and in radians per second. Returns graphics.";
PhasePlot::usage =
"PhasePlot[Transf, s, {w,wmin,wmax}, gopts]: Plots the phase part of Bode plot of transfer function Transf(s) in degree, from frequency w=wmin to w=wmax; w is in radian per second.Returns graphics.";
LinearLogPlot::usage =
"LinearLogPlot[f, {x, xmin, xmax}] generates a plot of f as a function of x \
in log scale.";
NyquistPlot::usage =
"NyquistPlot[Transf, s, {w,wmin,wmax}, gopts]: Plots the Nyquist plot of the transfer function Transf(s). The plot is composed of three parts. The first part corresponds to s=jw with w=-wmin to w=-wmax. The second part corresponds to s=jw, w varying from w=+wmin to w=+wmax. Encirclement information of (-1+j0) is provided by the third part which corresponds to s with theta=-pi/2 to pi/2. wmin and wmax should be > 0 and radians per second. Returns graphics.";
Polar::usage =
"Polar[Transf, s, {w,wmin,wmax}, gopts]: Plots the polar plot of the transfer function Transf(s) with s=jw , w varying from w=wmin to w=wmax, both > 0 and in radians per second. Returns graphics.";
MagvsPhase::usage =
"MagvsPhase[Transf, s, {w,wmin,wmax}, gopts]: Plots the magnitude vs. phase plot of transfer function Transf(s) with s=jw. The frequency w varies from w=wmin to w=wmax both > 0 and in radians per second. Returns graphics.";
Modify::usage=
"Modify[x] returns a proper value for phase angle";
Begin["`Private`"]
(* Magnitude plot of Bode plot *)
SetAttributes[MagPlot, HoldAll]
MagPlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{freqresp = Transf/. s->I w},
LinearLogPlot[20*Log[10,Abs[freqresp]],{w,wmin,wmax},Frame->True,
GridLines->{LogGridMinor,Automatic},
FrameLabel->{"Frequency (rad/sec)","20 log|G(s)| (dB)",
"Magnitude Plot"," "},opts]
]/; wmin < wmax
(* Modify angle of transfer function *)
Modify/:Modify[x_]:= x /; x < 0
Modify/:Modify[x_]:= -360 + x /; x > 0
(* Phase plot of Bode plot *)
SetAttributes[PhasePlot, HoldAll]
PhasePlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{freqresp = Transf/. s->I w},
LinearLogPlot[Modify[180*Arg[freqresp]/Pi],{w,wmin,wmax},
Frame->True,GridLines->{LogGridMinor,Automatic},
FrameLabel->{"Frequency (rad/sec)","Degree of |G(s)|",
"Phase Plot"," "},opts]
]/; wmin < wmax
(* dB vs. Phase Plot *)
MagvsPhase[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{freqresp = Transf/. s->I w},
ParametricPlot[Evaluate[{Modify[180*Arg[freqresp]/Pi],
20*Log[10,Abs[freqresp]]}],
{w, wmin, wmax},Ticks->{Automatic, Automatic},
Frame->True,GridLines->{Automatic,Automatic},
FrameTicks->{Automatic,Automatic,Automatic,Automatic},
FrameLabel->{"Degree of |G(s)|","20 log|G(s)| (dB)",
"dB vs. Phase Plot"," "},opts]
]/; wmin < wmax
(* Linear vs. Log plot *)
SetAttributes[LinearLogPlot, HoldAll]
LinearLogPlot[f_, {x_, xmin_, xmax_}, opts___] :=
Module[{},
ParametricPlot[{Log[10,x], f}, {x, xmin, xmax},
Ticks->{LogScale, Automatic},
FrameTicks -> {LogScale, Automatic,
LogScale, Automatic},
PlotRange -> All, opts]
]
(* Nyquist Plot *)
NyquistPlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{const,mag,phase,trans,npoint=40,gtable,theta,rad,theta1,
ang1,fac,wwmin,centerx,angle,firstpt,finalpt,p1,p2,
rad1,rad2,arg1,arg2,t,tmin,tmax,dfirst,dfinal,
gp,g,g1,gminus,glist,freqresp,num,den,nn,x,radius},
(* The frequency response *)
freqresp = Expand[Transf/. s->I w] ;
mag = Abs[freqresp];
phase = Arg[freqresp];
(* Plot the (-1,0) point *)
gp = Graphics[{PointSize[0.02],Point[{-1,0}]}];
glist={gp};
(* Manipulate the transfer function to standard form *)
trans = Together[ExpandAll[Transf]];
num = Numerator[trans];
den = Denominator[trans];
(* If there exists poles at zero *)
If[ Limit[den,s->0]==0,
(* Reduce the poles at s=0 *)
For[nn=0, Limit[den,s->0]==0, nn++, den=Simplify[den/s]];
wwmin = wmin;
const = (num/den/.s->0);
centerx = 0;
(* The type 1 problem *)
If[nn ==1 ,
firstpt = {mag Cos[phase] ,mag Sin[phase]}/.w->-wwmin;
finalpt = {mag Cos[phase],mag Sin[phase]}/.w->wwmin;
rad = Max[Abs[mag Cos[phase]],Abs[mag Sin[phase]]]/.w->-wwmin;
fac=1;
If[const<0, centerx = mag Cos[phase]/.w->-wwmin;
angle[x_]:= - Pi - x , angle[x_]:= - x];
theta[x_]:=-Pi/2 + Pi/npoint x;
gtable=Table[ ang1=angle[theta[i]];
{rad*Cos[ang1]+centerx, rad*Sin[ang1]},{i,0,npoint}];
gtable = Prepend[gtable,firstpt];
gtable = Append[gtable,finalpt];
(* Plot the graph from -0 to +0 theta from -90 to +90 *)
g1 = ListPlot[gtable, DisplayFunction -> Identity,
PlotJoined->True,
PlotStyle->{Dashing[{0.02,0.02}]}];
];
(* The type 2 & above problem *)
If[nn >= 2,
(* Factor to reduce the wwnin (to increase rad2) *)
fac = 0.95;
firstpt = {mag Cos[phase] ,mag Sin[phase]}/.w->-wwmin;
finalpt = {mag Cos[phase],mag Sin[phase]}/.w->wwmin*fac;
p1 = firstpt[[1]]+I firstpt[[2]];
p2 = finalpt[[1]]+I finalpt[[2]];
rad1 = Abs[ p1 ]; (* Initial radius *)
rad2 = Abs[ p2 ]; (* Final radius *)
(* Normalized initial angle *)
arg1=N[Mod[ Arg[ p1 ] ,Pi/2]];
(* Normalized final angle into 0<= arg <= Pi/2 *)
arg2=N[Mod[ Arg[ p2 ] ,Pi/2]];
If[arg1 < N[Pi/4], dfirst = arg1, dfirst = -Pi/2 + arg1];
If[arg2 < N[Pi/4], dfinal = -arg2, dfinal = Pi/2 - arg2];
If[const<0,
tmin= -Pi - nn(-Pi/2)+dfirst; tmax=-Pi - nn(Pi/2)-dfinal ,
tmin= - nn(-Pi/2)+dfirst; tmax= - nn(Pi/2)-dfinal];
radius[x_]:= rad1 + (x-tmin)/(tmax-tmin) (rad2 - rad1);
(* Plot the graph from -0 to +0 theta from -90 to +90 *)
g1 = PolarPlot[radius[t], {t,tmin,tmax},
DisplayFunction -> Identity,
PlotStyle->{Dashing[{0.02,0.02}]}]
];
glist = {gp,g1} ,
wwmin = 0 (* If there is no pole at zero, set wwmin = 0 *)
];
(* Plot the polar plot for w from -0 -> -inifinity and 0 ->+inifinity *)
g = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, {w, wwmin*fac, wmax},
DisplayFunction -> Identity, opts];
gminus = ParametricPlot[{mag Cos[phase], mag Sin[phase]},
{w, -wwmin , -wmax},DisplayFunction -> Identity];
glist = Join[glist,{g,gminus}];
Show[glist, DisplayFunction -> $DisplayFunction,
Ticks->{Automatic, Automatic},Frame->True,
FrameTicks -> {Automatic, Automatic,Automatic, Automatic},
FrameLabel->{"Re(G(jw))","Im(G(jw))","Nyquist Plot"," "},
GridLines->{Automatic,Automatic},
PlotRange -> All,opts
]
]/; wmin < wmax
(* Polar Plot *)
Polar[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
Module[{mag,phase,gp,g,glist,freqresp},
(* The frequency response *)
freqresp = Expand[Transf/. s->I w] ;
mag = Abs[freqresp];
phase = Arg[freqresp];
(* Plot the (-1,0) point *)
gp = Graphics[{PointSize[0.02],Point[{-1,0}]}];
(* Plot the polar plot for w from -0 -> -inifinity and 0 ->+inifinity *)
g = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, {w, wmin, wmax},
DisplayFunction -> Identity, opts];
glist={gp,g};
Show[glist, DisplayFunction -> $DisplayFunction,
Ticks->{Automatic, Automatic},Frame->True,
FrameTicks -> {Automatic, Automatic,Automatic, Automatic},
FrameLabel->{"Re(G(jw))","Im(G(jw))","Polar Plot"," "},
GridLines->{Automatic,Automatic},
PlotRange -> All,opts
]
]/; wmin < wmax
End[];
EndPackage[];
(* E N D P A C K A G E *)
If [ TrueQ[ $VersionNumber >= 2.0 ],
On[ General::spell ];
On[ General::spell1 ] ]
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.