--- a/main/control/inst/@ss/__c2d__.m	Fri Oct 26 17:07:34 2012 +0000
+++ b/main/control/inst/@ss/__c2d__.m	Sat Oct 27 11:58:21 2012 +0000
@@ -45,6 +45,14 @@
         [sys.a, sys.b, sys.c, sys.d, sys.e] = __dss_bilin__ (sys.a, sys.b, sys.c, sys.d, sys.e, beta, false);
       endif
 
+    case "m"                       # "matched"
+      tmp = ss (c2d (zpk (sys), tsam, method));
+      sys.e = tmp.e;
+      sys.a = tmp.a;
+      sys.b = tmp.b;
+      sys.c = tmp.c;
+      sys.d = tmp.d;
+
     otherwise
       error ("ss: c2d: %s is an invalid or missing method", method);
   endswitch

--- a/main/control/inst/@tf/__c2d__.m	Fri Oct 26 17:07:34 2012 +0000
+++ b/main/control/inst/@tf/__c2d__.m	Sat Oct 27 11:58:21 2012 +0000
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2011   Lukas F. Reichlin
+## Copyright (C) 2009, 2011, 2012   Lukas F. Reichlin
 ##
 ## This file is part of LTI Syncope.
 ##
@@ -20,22 +20,60 @@
 
 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 ## Created: October 2009
-## Version: 0.2
+## Version: 0.3
 
 function sys = __c2d__ (sys, tsam, method = "zoh", w0 = 0)
 
-  [p, m] = size (sys);
+  if (strncmpi (method, "m", 1))    # "matched"
+    ## TODO: move this code to  @zpk/__c2d__.m  once ZPK models are implemented
+    
+    if (! issiso (sys))
+      error ("tf: c2d: require SISO system for matched pole/zero method");
+    endif
+
+    [z_c, p_c, k_c] = zpkdata (sys, "vector");
+    p_d = exp (p_c * tsam);    
+    z_d = exp (z_c * tsam);
+
+    if (any (! isfinite (p_d)) || any (! isfinite (z_d)))
+      error ("tf: c2d: discrete-time poles and zeros are not finite");
+    endif
+
+    ## continuous-time zeros at infinity are mapped to -1 in discrete-time
+    ## except for one.  for non-proper transfer functions, no zeros at -1 are added.
+    np = length (p_c);              # number of poles
+    nz = length (z_c);              # number of finite zeros, np-nz number of infinite zeros
+    z_d = vertcat (z_d, repmat (-1, np-nz-1, 1));
 
-  for i = 1 : p
-    for j = 1 : m
-      idx = substruct ("()", {i, j});
-      tmp = subsref (sys, idx);
-      tmp = c2d (ss (tmp), tsam, method, w0);
-      [num, den] = tfdata (tmp, "tfpoly");
-      sys.num(i, j) = num;
-      sys.den(i, j) = den;
+    ## the discrete-time gain k_d is matched at a certain frequency (w_c, w_d)
+    ## to continuous-time gain k_c.  the dc gain is taken (w_c=0, w_d=1) unless
+    ## there are continuous-time poles/zeros near 0.  then w_c=1/tsam is taken.
+    w_c = 0;                        # dc gain
+    tol = sqrt (eps);               # poles/zeros below tol are assumed to be zero
+    while (any (abs ([p_c; z_c] - w_c) < tol))
+      w_c += 0.1 / tsam;
+    endwhile
+    w_d = exp (w_c * tsam);
+    k_d = real (k_c * prod (w_c - z_c) / prod (w_c - p_c) * prod (w_d - p_d) / prod (w_d - z_d));
+
+    tmp = zpk (z_d, p_d, k_d, tsam);
+    sys.num = tmp.num;
+    sys.den = tmp.den;
+
+  else
+    [p, m] = size (sys);
+
+    for i = 1 : p
+      for j = 1 : m
+        idx = substruct ("()", {i, j});
+        tmp = subsref (sys, idx);
+        tmp = c2d (ss (tmp), tsam, method, w0);
+        [num, den] = tfdata (tmp, "tfpoly");
+        sys.num(i, j) = num;
+        sys.den(i, j) = den;
+      endfor
     endfor
-  endfor
+  endif
 
   sys.tfvar = "z";