getf SCI/util/testexamples.sci
reinit_for_test()
%U=mopen('SCI/tests/automatic_tests/sfact_data.ref','rb');
//Simple polynomial example
z = poly(0, 'z');
p = (z - 1/2) * (2 - z);
if load_ref('p') then   pause,end,

w = sfact(p);
%ans = w * numer(horner(w, 1/z));
if load_ref('%ans') then   pause,end,

//matrix example
F1 = [z - 1/2,z + 1/2,z^2 + 2;1,z,-z;z^3 + 2 * z,z,1/2 - z];
P = F1 * gtild(F1, 'd');//P is symmetric
F = sfact(P);
if load_ref('F') then   pause,end,

%ans = roots(det(P));
if load_ref('%ans') then   pause,end,

%ans = roots(det(gtild(F, 'd')));
if load_ref('%ans') then   pause,end,
//The stable roots
%ans = roots(det(F));
if load_ref('%ans') then   pause,end,
//The antistable roots
%ans = clean(P - F * gtild(F, 'd'));
if load_ref('%ans') then   pause,end,

//Example of continuous time use
s = poly(0, 's');
p = -3 * (s + 1 + %i) * (s + 1 - %i) * (s + 0.5) * (s - 0.5) * (s - (1 + %i)) * (s - (1 - %i));p = real(p);
//p(s) = polynomial in s^2 , looks for stable f such that p=f(s)*f(-s)
w = horner(p, (1 - s)/(1 + s));// bilinear transform w=p((1-s)/(1+s))
wn = numer(w);//take the numerator
fn = sfact(wn);f = numer(horner(fn, (1 - s)/(s + 1)));//Factor and back transform
f = f/sqrt(horner(f * gtild(f, 'c'), 0));f = f * sqrt(horner(p, 0));//normalization
%ans = roots(f);
if load_ref('%ans') then   pause,end,
//f is stable
%ans = clean(f * gtild(f, 'c') - p);
if load_ref('%ans') then   pause,end,
//f(s)*f(-s) is p(s)
xdel_run(winsid());

mclose(%U);