# HG changeset patch
# User octavecontrib@gmail.com
# Date 1536074609 -19800
#      Tue Sep 04 20:53:29 2018 +0530
# Node ID 90dbd48910b15116e2b2e95d996c7896c64e3437
# Parent  fa91b35bee7151b8d81cfcb92d6d0ce6454f98d8
Improvement to frequency point selection for bode, nyquist, nichols etc. Small change to makefile.

* __frequency_vector__.m: A new local function '(__w_based_on_complex_conj_roots__)'
produces a vector of frequency points. This is then passed to '(__frequency_vector__)'
function through '(__frequency_range__)' function. This vector is merged into the
existing logspaced frequency points. Some tests are added to the end of this file.
* Makefile
Writes the existing PKG_ADD to PKG_ADD file in src directory.
Removed --no-window-system option from the default list.

diff -r fa91b35bee71 -r 90dbd48910b1 Makefile
--- a/Makefile	Fri May 04 04:10:35 2018 +0200
+++ b/Makefile	Tue Sep 04 20:53:29 2018 +0530
@@ -19,7 +19,7 @@ M_SOURCES   := $(wildcard inst/*.m)
 CC_SOURCES  := $(wildcard src/*.cc)
 PKG_ADD     := $(shell grep -sPho '(?<=(//|\#\#) PKG_ADD: ).*' $(CC_SOURCES) $(M_SOURCES))
 
-OCTAVE ?= octave --no-window-system --silent
+OCTAVE ?= octave --silent
 
 .PHONY: help dist html release install all check run clean
 
@@ -70,19 +70,19 @@ release: dist html
 ## dependencies on DESCRIPTION.
 install: $(RELEASE_TARBALL)
 	@echo "Installing package locally ..."
-	$(OCTAVE) --eval 'pkg ("install", "${RELEASE_TARBALL}")'
+	$(OCTAVE) --no-window-system --eval 'pkg ("install", "${RELEASE_TARBALL}")'
 
 all: $(CC_SOURCES)
 	$(MAKE) -C src/ all
+	@echo '${PKG_ADD}' > src/PKG_ADD
 
 check: all
-	$(OCTAVE) --path "inst/" --path "src/" \
+	$(OCTAVE) --no-window-system --path "inst/" --path "src/" \
 	  --eval '${PKG_ADD}' \
 	  --eval 'runtests ("inst"); runtests ("src");'
 
 run: all
-	$(OCTAVE) --persist --path "inst/" --path "src/" \
-	  --eval '${PKG_ADD}'
+	$(OCTAVE) --persist --path "inst/" --path "src/"
 
 clean:
 	$(RM) -r $(TARGET_DIR)
diff -r fa91b35bee71 -r 90dbd48910b1 inst/__frequency_vector__.m
--- a/inst/__frequency_vector__.m	Fri May 04 04:10:35 2018 +0200
+++ b/inst/__frequency_vector__.m	Tue Sep 04 20:53:29 2018 +0530
@@ -49,9 +49,12 @@ function w = __frequency_vector__ (sys_c
   sys_cell = sys_cell(idx);
   len = numel (sys_cell);
   
-  [dec_min, dec_max, zp] = cellfun (@__frequency_range__, sys_cell, {wbounds}, "uniformoutput", false);
+  [dec_min, dec_max, zp, addlw] = cellfun (@__frequency_range__, ...
+                                           sys_cell, {wbounds}, ...
+                                           "uniformoutput", false ...
+                                          );
 
-  if (strcmpi (wbounds, "std"))                 # plots with explicit frequencies
+  if (strcmpi (wbounds, "std"))             # plots with explicit frequencies
 
     if (nargin == 4)                        # w = {wmin, wmax}  
       dec_min = log10 (wmin);
@@ -67,8 +70,15 @@ function w = __frequency_vector__ (sys_c
     idx = find (zp > 10^dec_min & zp < 10^dec_max);
     zp = zp(idx);
 
-    w = logspace (dec_min, dec_max, 500);
-    w = unique ([w, zp]);                       # unique also sorts frequency vector
+    addlw = horzcat (addlw{:});
+
+    ## remove additional frequency points outside (dec_min, dec_max)
+    idx = find (addlw > 10^dec_min & addlw < 10^dec_max);
+    addlw = addlw(idx);
+
+
+    w = logspace (dec_min, dec_max, 50);
+    w = unique ([w, zp, addlw]);         # Unique also sorts frequency vector.
 
     w = repmat ({w}, 1, len);                   # return cell of frequency vectors
 
@@ -81,10 +91,18 @@ function w = __frequency_vector__ (sys_c
 
     idx = cellfun (@(zp, dec_min, dec_max) find (zp > 10^dec_min & zp < 10^dec_max), ...
                    zp, dec_min, dec_max, "uniformoutput", false);
-    zp = cellfun (@(zp, idx) zp(idx), zp, idx, "uniformoutput", false);
+    zp  = cellfun (@(zp, idx) zp(idx), zp, idx, "uniformoutput", false);
 
-    w = cellfun (@logspace, dec_min, dec_max, {500}, "uniformoutput", false);
-    w = cellfun (@(w, zp) unique ([w, zp]), w, zp, "uniformoutput", false);
+    idx = cellfun (@(addlw, dec_min, dec_max) find (addlw > 10^dec_min ...
+                                                    & addlw < 10^dec_max), ...
+                   addlw, dec_min, dec_max, "uniformoutput", false);
+    addlw = cellfun (@(addlw, idx) addlw(idx), ...
+                     addlw, idx, "uniformoutput", false);
+
+
+    w = cellfun (@logspace, dec_min, dec_max, {50}, "uniformoutput", false);
+    w = cellfun (@(w, zp, addlw) unique ([w, zp, addlw]), ...
+                 w, zp, addlw, "uniformoutput", false);
     ## unique also sorts frequency vector
 
   else
@@ -98,13 +116,14 @@ function w = __frequency_vector__ (sys_c
 endfunction
 
 
-function [dec_min, dec_max, zp] = __frequency_range__ (sys, wbounds = "std")
+function [dec_min, dec_max, zp, addlw] = __frequency_range__ (sys, wbounds = "std")
 
   if (isa (sys, "frd"))
     w = get (sys, "w");
     dec_min = log10 (w(1));
     dec_max = log10 (w(end));
     zp = [];
+    addlw = [];
     return;
   endif
     
@@ -191,4 +210,164 @@ function [dec_min, dec_max, zp] = __freq
   ## include zeros and poles for nice peaks in plots
   zp = [abs(zer), abs(pol)];
 
+  addlw = [];
+  if (! discrete)
+    addlw = __w_based_on_complex_conj_roots__ ([zer, pol]);
+  endif
+
 endfunction
+
+function w = __w_based_on_complex_conj_roots__ (zp)
+  ## Returns a list of frequency points based on
+  ## magnitude and phase contribution from
+  ## complex roots (poles and zeroes).
+  ## The formula is applicable for continuous time systems.
+  ## The formula doesn't account for multiplicity of roots.
+  ## The formula doesn't account for the residue of roots
+  ## in the partial fraction expansion of the
+  ## system transfer function.
+  ##
+  ## Input is expected to be a row vector of complex numbers.
+  ## Output is a row vector.
+
+  w = [];
+
+  if ((! isvector (zp)) || rows (zp) != 1 || length (zp) == 0)
+    return;
+  endif
+
+  ## Filter the input for one root from a complex conjugate pair.
+  ## Remove purely real roots and roots at origin.
+  zp = zp( imag (zp) > eps );
+
+  ## Remove roots at infinity.
+  zp = zp( abs (zp)  < Inf );
+
+
+  ## #######################################
+  ## Calculation based on phase contribution
+  ## #######################################
+
+  ph  = [-1,  -5 : -5 : -175, -179] *pi/180; # in radians, configurable.
+  tph = tan (ph);
+  ## Above ph and its tangent (tan (ph)) may be better off
+  ## being persistent variables.
+
+  bb = @(x) 2 * real (x) ./ tph;
+  cc = @(x) -x * x';
+
+  ## Solution for quadratic equation to be applied for
+  ## all roots supplied.
+  w1 = arrayfun (@(x)(-bb (x) + sqrt ((bb (x)).^2 - 4 * cc (x))) / 2, ...
+                 zp, "UniformOutput", false ...
+                );
+  w1 = horzcat (w1{:});
+
+  clear bb cc
+
+  ## ##############################
+  ## Calculation based on magnitude
+  ## ##############################
+
+  yy_at_zero = 1 ./ (abs (zp)) .^ 2;                 # mag at zero freq
+  yy_at_wn   = 1 ./ abs (2 * real (zp) .* abs (zp)); # mag at natural freq
+  yy_low     =  min (yy_at_zero, yy_at_wn);
+  yy_high    =  max (yy_at_zero, yy_at_wn);
+  decades    = ceil (log10 (yy_high ./ yy_low));
+  decades    =  min (decades, 50);                # remove Inf if any
+  yy_high    =  min (yy_high, yy_low .* 10^50);   # remove Inf if any
+  points_per_decade = 4;                          # configurable number
+
+  yy = arrayfun (@logspace, log10 (yy_low), log10 (yy_high), ...
+                 points_per_decade .* decades, ...
+                 "UniformOutput", false ...
+                );
+
+  bb = @(x) 4 * (real (x)) .^ 2 - 2 * (abs (x)) .^ 2;
+  cc = @(x, y) (abs (x)) .^ 4 - 1 ./ y .^ 2;
+
+  ## Solution for quadratic equation to be applied for
+  ## all roots supplied.
+  w2sqd = cellfun (@(x, y) (-bb (x) + sqrt ((bb (x)) .^ 2 ...
+                                            - 4 * cc (x, y))) / 2, ...
+                   num2cell (zp), yy, "UniformOutput", false ...
+                  );
+  w3sqd = cellfun (@(x, y) (-bb (x) - sqrt ((bb (x)) .^ 2 ...
+                                            - 4 * cc (x, y))) / 2, ...
+                   num2cell (zp), yy, "UniformOutput", false ...
+                  );
+
+  w2sqd = horzcat (w2sqd{:});
+  w3sqd = horzcat (w3sqd{:});
+
+  ## Since w2sqd and w3sqd give the square of frequency points
+  ## on either side of the frequency of the root,
+  ## a solution need not exist for roots with no 'M-peak' in
+  ## in their trasfer function.
+  ## So, only positve real numbers whose sqrt are positve real
+  ## need to be considered for further processing.
+  w2sqd = w2sqd( real (w2sqd) >= 0 & abs (imag (w2sqd)) < eps);
+  w3sqd = w3sqd( real (w3sqd) >= 0 & abs (imag (w3sqd)) < eps);
+
+  w2 = sqrt (w2sqd);
+  w3 = sqrt (w3sqd);
+
+  ## Do sanity check on the frequency points.
+  ## FIXME: To move these checks to test blocks.
+  assert (all (abs (imag (w1)) <  eps), "w1 has imaginary part.");
+  assert (all (real (w1)  >=   0),      "w1 has negative real part.");
+  assert (! any (isnan (w1)),           "w1 contains NaN.");
+
+  assert (all (abs (imag (w2)) <  eps), "w2 has imaginary part.");
+  assert (all (real (w2)  >=   0),      "w2 has negative real part.");
+  assert (! any (isnan (w2)),           "w2 contains NaN.");
+
+  assert (all (abs (imag (w3)) <  eps), "w3 has imaginary part.");
+  assert (all (real (w3)  >=   0),      "w3 has negative real part.");
+  assert (! any ( isnan (w3)),          "w3 contains NaN.");
+
+  ## Join all into one vector.
+  w  = unique ([w1, w2, w3]);
+endfunction
+
+%!test
+%! sys{1}  = tf (         [1],          [1]);
+%! sys{2}  = tf (         [1],      [1,  1]);
+%! sys{3}  = tf (     [1,  1],          [1]);
+%! sys{4}  = tf (     [1,  1],      [1,  3]);
+%! sys{5}  = tf (     [1,  3],      [1,  1]);
+%! sys{6}  = tf (         [1], [1, 0.1, 10]);
+%! sys{7}  = tf ([1,   1, 10],          [1]);
+%! sys{8}  = tf ([1,   3, 10], [1,   1, 10]);
+%! sys{9}  = tf ([1,   1, 10], [1,   3, 10]);
+%! sys{10} = tf ([1, 0.1,  1], [1,   1,  1]) * ...
+%!           tf ([1, 100, 1e8],[1, 300, 1e8]);
+%! sys{11} = tf ([1, 0.1,  1], [1,   1,  1]) * ...
+%!           tf ([1, 100, 1e8], [1, 300, 1e8]) * tf ([1], [1, 0]);
+%! sys{12} = tf (         [1],      [1,  0]);
+%! sys{13} = tf (     [1,  0],      [1,100]);
+%! sys{14} = sys{6} * sys{6};
+%! sys{15} = sys{6} * sys{6} * sys{6};
+%! sys{16} = tf ([1,  0.01, 10], [1, 0.01, 15]);
+%! sys{17} = sys{16} * tf ([1,  0.01, 20], [1, 0.01, 25]);
+%! sys{18} = tf ([1,  0.001, 150000], [1, 0.001, 160000]);
+%! sys{19} = sys{18} * tf ([1,  0.001, 170000], [1, 0.001, 180000]);
+%! sys{20} = sys{17} * sys{19};
+%! #
+%! # Poles on imaginary axis
+%! sys{21} = tf ([1], [1, 0, 1]);
+%! sys{22} = tf ([1, 0, 100], [1, 0, 10000]);
+%! #
+%! # Roots on right side
+%! sys{23} = tf ([1,  -1, 10],          [1]);
+%! sys{24} = tf ([1,  -3, 10], [1,  -1, 10]);
+%! sys{25} = tf ([1,  -1, 10], [1,  -3, 10]);
+%! # All pass filter
+%! sys{26} = tf ([-1, 1], [1, 1]);
+%! #
+%! # Start actual testing
+%! [w] = cellfun (@__frequency_vector__, sys, "uniformoutput", false);
+%! [w] = __frequency_vector__ (sys, "std");
+%! [w] = __frequency_vector__ (sys, "ext");
+%! [w] = __frequency_vector__ (sys, "std", {1, 50});
+%! [w] = __frequency_vector__ (sys, "ext", {1, 50});
