home *** CD-ROM | disk | FTP | other *** search
/ Education Sampler 1992 [NeXTSTEP] / Education_1992_Sampler.iso / Mathematics / Notebooks / COSY_PAK / COSYPAK / chap6.m < prev    next >
Text File  |  1992-06-02  |  9KB  |  234 lines
  1.  
  2. (* 
  3. **********************************************************
  4.  
  5.          COSY_PAK package chap6.m        
  6.  
  7. **********************************************************
  8. *)
  9.  
  10.  
  11. If [ TrueQ[ $VersionNumber >= 2.0 ],
  12.      Off[ General::spell ];
  13.      Off[ General::spell1 ] ]
  14.      
  15. (*  B E G I N     P A C K A G E  *)
  16.  
  17.  
  18. BeginPackage["COSYPAK`chap6`","Algebra`ReIm`",
  19.          "Graphics`Graphics`"]
  20.          
  21.          
  22. MagPlot::usage = 
  23.  "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.";  
  24.            
  25. PhasePlot::usage =
  26. "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.";  
  27.  
  28. LinearLogPlot::usage =
  29. "LinearLogPlot[f, {x, xmin, xmax}] generates a plot of f as a function of x \
  30.  in log scale.";
  31.  
  32. NyquistPlot::usage =
  33. "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.";
  34.  
  35. Polar::usage =
  36. "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.";
  37.  
  38. MagvsPhase::usage = 
  39. "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.";
  40.  
  41. Modify::usage=
  42. "Modify[x] returns a proper value for phase angle";
  43.  
  44. Begin["`Private`"]
  45.  
  46. (* Magnitude plot of Bode plot *)
  47. SetAttributes[MagPlot, HoldAll]
  48.  
  49. MagPlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=         
  50. Module[{freqresp = Transf/. s->I w},
  51.     LinearLogPlot[20*Log[10,Abs[freqresp]],{w,wmin,wmax},Frame->True,
  52.     GridLines->{LogGridMinor,Automatic},
  53.     FrameLabel->{"Frequency (rad/sec)","20 log|G(s)| (dB)",
  54.     "Magnitude Plot"," "},opts]
  55.     ]/; wmin < wmax    
  56.  
  57. (* Modify angle of transfer function *)
  58. Modify/:Modify[x_]:= x /; x < 0
  59. Modify/:Modify[x_]:= -360 + x /; x > 0
  60.  
  61. (* Phase plot of Bode plot *)
  62. SetAttributes[PhasePlot, HoldAll]
  63.  
  64. PhasePlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
  65. Module[{freqresp = Transf/. s->I w},
  66.     LinearLogPlot[Modify[180*Arg[freqresp]/Pi],{w,wmin,wmax},
  67.         Frame->True,GridLines->{LogGridMinor,Automatic},
  68.         FrameLabel->{"Frequency (rad/sec)","Degree of |G(s)|",
  69.         "Phase Plot"," "},opts]
  70.     ]/; wmin < wmax    
  71.  
  72. (* dB vs. Phase Plot *)
  73. MagvsPhase[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
  74. Module[{freqresp = Transf/. s->I w},
  75.      ParametricPlot[Evaluate[{Modify[180*Arg[freqresp]/Pi],
  76.              20*Log[10,Abs[freqresp]]}], 
  77.             {w, wmin, wmax},Ticks->{Automatic, Automatic},
  78.             Frame->True,GridLines->{Automatic,Automatic}, 
  79.             FrameTicks->{Automatic,Automatic,Automatic,Automatic},
  80.             FrameLabel->{"Degree of |G(s)|","20 log|G(s)| (dB)",
  81.             "dB  vs. Phase Plot"," "},opts]       
  82.       ]/; wmin < wmax
  83.  
  84. (* Linear vs. Log plot *)
  85. SetAttributes[LinearLogPlot, HoldAll]
  86.  
  87. LinearLogPlot[f_, {x_, xmin_, xmax_}, opts___] :=
  88.     Module[{},
  89.          ParametricPlot[{Log[10,x], f}, {x, xmin, xmax},
  90.                      Ticks->{LogScale, Automatic}, 
  91.                      FrameTicks -> {LogScale, Automatic,
  92.                            LogScale, Automatic},
  93.                     PlotRange -> All, opts]
  94.            ]
  95.        
  96.       
  97. (* Nyquist Plot *)    
  98. NyquistPlot[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
  99. Module[{const,mag,phase,trans,npoint=40,gtable,theta,rad,theta1,
  100.         ang1,fac,wwmin,centerx,angle,firstpt,finalpt,p1,p2,
  101.         rad1,rad2,arg1,arg2,t,tmin,tmax,dfirst,dfinal,
  102.         gp,g,g1,gminus,glist,freqresp,num,den,nn,x,radius},    
  103.     (* The frequency response *)
  104.         freqresp = Expand[Transf/. s->I w] ;
  105.             mag = Abs[freqresp]; 
  106.         phase = Arg[freqresp];
  107.     (* Plot the (-1,0) point *)
  108.         gp = Graphics[{PointSize[0.02],Point[{-1,0}]}];
  109.         glist={gp};
  110.     (* Manipulate the transfer function to standard form *)    
  111.         trans = Together[ExpandAll[Transf]];    
  112.         num = Numerator[trans];
  113.         den = Denominator[trans];
  114.     (* If there exists poles at zero *)        
  115.         If[ Limit[den,s->0]==0,
  116.            (* Reduce the poles at s=0 *)
  117.             For[nn=0, Limit[den,s->0]==0, nn++, den=Simplify[den/s]];
  118.             wwmin = wmin;
  119.             const = (num/den/.s->0);
  120.               centerx = 0;
  121.             
  122.             
  123.             
  124.         (* The type 1 problem *)
  125.             If[nn ==1 , 
  126.               firstpt = {mag Cos[phase] ,mag Sin[phase]}/.w->-wwmin;        
  127.               finalpt = {mag Cos[phase],mag Sin[phase]}/.w->wwmin;
  128.  
  129.               rad = Max[Abs[mag Cos[phase]],Abs[mag Sin[phase]]]/.w->-wwmin;
  130.               fac=1;              
  131.               If[const<0, centerx = mag Cos[phase]/.w->-wwmin;
  132.                        angle[x_]:= - Pi -  x , angle[x_]:= - x];
  133.               theta[x_]:=-Pi/2 + Pi/npoint x;            
  134.               gtable=Table[ ang1=angle[theta[i]];
  135.                 {rad*Cos[ang1]+centerx, rad*Sin[ang1]},{i,0,npoint}];
  136.             
  137.               gtable = Prepend[gtable,firstpt];
  138.               gtable = Append[gtable,finalpt];    
  139.             (* Plot the graph from -0 to +0 theta from -90 to +90 *)    
  140.               g1 = ListPlot[gtable, DisplayFunction -> Identity,
  141.                 PlotJoined->True,
  142.                 PlotStyle->{Dashing[{0.02,0.02}]}];
  143.          
  144.                        
  145.               ];
  146.               
  147.         (* The type 2 & above problem *)  
  148.             If[nn >= 2, 
  149.                (* Factor to reduce the wwnin (to increase rad2) *)
  150.                fac = 0.95; 
  151.                firstpt = {mag Cos[phase] ,mag Sin[phase]}/.w->-wwmin;        
  152.                finalpt = {mag Cos[phase],mag Sin[phase]}/.w->wwmin*fac;
  153.                p1 = firstpt[[1]]+I firstpt[[2]];
  154.                p2 = finalpt[[1]]+I finalpt[[2]];
  155.                rad1 = Abs[ p1 ]; (* Initial radius *)
  156.                rad2 = Abs[ p2 ]; (* Final radius *)
  157.                       (* Normalized initial angle *)
  158.                arg1=N[Mod[ Arg[ p1 ] ,Pi/2]]; 
  159.                (* Normalized final angle into 0<= arg <= Pi/2 *)
  160.                arg2=N[Mod[ Arg[ p2 ] ,Pi/2]];                
  161.                If[arg1 < N[Pi/4], dfirst = arg1, dfirst = -Pi/2 + arg1]; 
  162.                If[arg2 < N[Pi/4], dfinal = -arg2, dfinal = Pi/2 - arg2];
  163.        
  164.                If[const<0,
  165.                 tmin= -Pi - nn(-Pi/2)+dfirst; tmax=-Pi - nn(Pi/2)-dfinal , 
  166.                 tmin= - nn(-Pi/2)+dfirst; tmax= - nn(Pi/2)-dfinal];
  167.                            
  168.                radius[x_]:= rad1 + (x-tmin)/(tmax-tmin) (rad2 - rad1);
  169.                (* Plot the graph from -0 to +0 theta from -90 to +90 *)
  170.                g1 = PolarPlot[radius[t], {t,tmin,tmax}, 
  171.                DisplayFunction -> Identity,
  172.                PlotStyle->{Dashing[{0.02,0.02}]}]
  173.              ];   
  174.         
  175.             glist = {gp,g1} ,
  176.             
  177.             wwmin = 0     (* If there is no pole at zero, set wwmin = 0 *)    
  178.           ];
  179.  
  180.     (* Plot the polar plot for w from -0 -> -inifinity  and  0 ->+inifinity  *)
  181.         g = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, {w, wwmin*fac, wmax},
  182.                  DisplayFunction -> Identity, opts];
  183.         
  184.         gminus = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, 
  185.          {w, -wwmin  , -wmax},DisplayFunction -> Identity];
  186.         glist = Join[glist,{g,gminus}];
  187.         
  188.  
  189.         Show[glist, DisplayFunction -> $DisplayFunction,
  190.             Ticks->{Automatic, Automatic},Frame->True, 
  191.                 FrameTicks -> {Automatic, Automatic,Automatic, Automatic},
  192.         FrameLabel->{"Re(G(jw))","Im(G(jw))","Nyquist Plot"," "},
  193.         GridLines->{Automatic,Automatic}, 
  194.         PlotRange -> All,opts 
  195.         ]
  196.     ]/; wmin < wmax
  197.  
  198. (* Polar Plot *)    
  199. Polar[Transf_,s_,{w_,wmin_,wmax_},opts___] :=
  200. Module[{mag,phase,gp,g,glist,freqresp},    
  201.     (* The frequency response *)
  202.         freqresp = Expand[Transf/. s->I w] ;
  203.             mag = Abs[freqresp]; 
  204.         phase = Arg[freqresp];
  205.  
  206.     (* Plot the (-1,0) point *)
  207.         gp = Graphics[{PointSize[0.02],Point[{-1,0}]}];
  208.  
  209.     (* Plot the polar plot for w from -0 -> -inifinity  and  0 ->+inifinity  *)
  210.         g = ParametricPlot[{mag Cos[phase], mag Sin[phase]}, {w, wmin, wmax},
  211.                  DisplayFunction -> Identity, opts];
  212.  
  213.     glist={gp,g};         
  214.         
  215.     Show[glist, DisplayFunction -> $DisplayFunction,
  216.             Ticks->{Automatic, Automatic},Frame->True, 
  217.                 FrameTicks -> {Automatic, Automatic,Automatic, Automatic},
  218.         FrameLabel->{"Re(G(jw))","Im(G(jw))","Polar Plot"," "},
  219.         GridLines->{Automatic,Automatic}, 
  220.         PlotRange -> All,opts 
  221.         ]
  222.  
  223.     ]/; wmin < wmax
  224.  
  225.  
  226. End[];    
  227. EndPackage[];                  
  228.       
  229. (*  E N D     P A C K A G E  *)
  230.  
  231. If [ TrueQ[ $VersionNumber >= 2.0 ],
  232.      On[ General::spell ];
  233.      On[ General::spell1 ] ]
  234.