home *** CD-ROM | disk | FTP | other *** search
/ Education Sampler 1992 [NeXTSTEP] / Education_1992_Sampler.iso / Mathematics / Notebooks / COSY_PAK / COSYPAK / chap5.m < prev    next >
Text File  |  1992-05-31  |  4KB  |  121 lines
  1.  
  2. (* 
  3. **********************************************************
  4.  
  5.          COSY_PAK package chap5.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`chap5`",
  19.         "Graphics`Graphics`"];
  20.  
  21. RootLocus::usage = 
  22. "RootLocus[Transf, s, {k,kmin,kmax}, gopts]: Plots the root-locus plot of the transfer function Transf(s) with the Parameter k varying from k=kmin to k=kmax. ";
  23.  
  24. Begin["`Private`"];
  25.  
  26. (* Plot the Root Locus of a transfer function *)
  27. RootLocus[transf_,s_,{Para_,Paramin_,Paramax_},opts___]:=
  28. Module[{CharPoly,denom,numer,k,klist,poles,npoles,zeros,start,
  29.         g,glist,realline,grealline,sol,linelem,linelist,nlines},
  30.     
  31.     denom = Denominator[transf];
  32.     numer = Numerator[transf];
  33.     
  34.     (* Get the starting pole location *)
  35.     poles=NSolve[denom==0,s];
  36.     realline = Cases[s/.poles, x_/;(Abs[Im[x]] < 0.0001) ]; 
  37.     poles=Map[{Re[#],Im[#]}& ,s/.poles];
  38.         
  39.     poles=Table[Append[poles[[i]],"X"],{i,1,Length[poles]}];
  40.     npoles=Length[poles];
  41.     glist = Array[g,npoles+1,0];
  42.     start = poles;
  43.  
  44.     
  45.     (* If numerator exists *)    
  46.     If[! FreeQ[numer,s],    
  47.         
  48.         (* Produce the mark of zeros *)
  49.         zeros=NSolve[numer==0,s];
  50.         realline = Join[realline,Cases[s/.zeros, x_/;(Abs[Im[x]]<0.0001)]];  
  51.         zeros=Map[{Re[#],Im[#]}& ,s/.zeros];
  52.         zeros=Table[Append[zeros[[i]],"O"],{i,1,Length[zeros]}];
  53.         start=Join[poles,zeros]
  54.         ];
  55.     
  56.     (* Procedure to produce lines to push root locus to reach zero *) 
  57.     nlines = Ceiling[Length[realline]/2];
  58.        realline = Sort[realline,Greater];
  59.         linelist = Array[linelem, nlines ];
  60.     
  61.         If[EvenQ[Length[realline]],     
  62.              Do[ linelem[i] = Line[ {{realline[[2 i-1]],0},{realline[[2 i]],0}}],
  63.             {i,1, nlines}],
  64.             
  65.               Do[ linelem[i] = Line[ {{realline[[2 i-1]],0},{realline[[2 i]],0}} ],
  66.                 {i,1, nlines-1}];
  67.       
  68.           linelem[nlines] = 
  69.           Line[{{realline[[2 nlines-1]],0},{realline[[2 nlines -1]] -2 ,0}}] 
  70.            
  71.         ];
  72.       
  73.     linelist = Prepend[linelist,Thickness[0.0051]];
  74.     grealline = Graphics[linelist]; 
  75.     CharPoly = Collect[Numerator[Together[1 + transf Para]],s];
  76.     
  77.     (* Looking for positive break away or break in gain *)
  78.     keqn = Para/.Solve[CharPoly == 0,Para][[1]];
  79.     sol = NSolve[Numerator[Together[D[keqn,s]]]==0,s];
  80.  
  81.     (* Pick up positive real parameter at the breaking point *)
  82.     klist = Cases[keqn/.sol,x_/;(Re[x]>0 && Abs[Im[x]] < 0.0001) ];
  83.  
  84.     (* Producing the parameter listing for plotting *)
  85.     k = Paramin;
  86.     a = 0.08;b = 0.03;
  87.     While[k <= Paramax ,klist = Append[klist,k]; k = k (1 + a) + b];
  88.  
  89.     (* Appedn the final value of k then sort the klist *)
  90.     klist = Sort[Append[klist,Paramax]];
  91.     EqnSolve[x_]:= NSolve[ (CharPoly/.Para->x) == 0 ,s];
  92.     points = Map[{Re[#],Im[#]}& , s/.Map[EqnSolve,klist] ,{2}]; 
  93.     (* Graph of Poles and Zeros position *)
  94.     g[0] = TextListPlot[start,DisplayFunction -> Identity];
  95.  
  96.  
  97.     (* Graph Sections of Root Locus devided by break points *)
  98.     Do[ g[j] = ListPlot[ Table[points[[i,j]],{i,1,Length[klist]}],
  99.             PlotJoined->True,PlotStyle->{Thickness[0.0051]},
  100.             DisplayFunction -> Identity ],
  101.         {j,1,npoles}
  102.       ];
  103.     (* Show all graphics *)    
  104.     Show[{grealline,glist},
  105.         DisplayFunction -> $DisplayFunction,PlotRange -> All,
  106.         FrameLabel->{"Re","Im","Root Locus"," "},Frame->True,
  107.         Ticks->{Automatic, Automatic},
  108.         GridLines->Automatic, 
  109.         FrameTicks -> {Automatic, Automatic,Automatic, Automatic},opts]
  110.     ]/; Paramin < Paramax
  111.         
  112. End[];
  113. EndPackage[];
  114.  
  115. (*  E N D     P A C K A G E  *)
  116.  
  117. If [ TrueQ[ $VersionNumber >= 2.0 ],
  118.      On[ General::spell ];
  119.      On[ General::spell1 ] ]
  120.  
  121.