% Discretization using Tustin's scheme via CFE (continued fractional expansion)
% using MAPLE under MATLAB6. 26/12/00
% by Dr YangQuan Chen <yqchen@ieee.org>
% function [b,a]=tustin_cfe(n,r,ts)
% b: numerator; a: denominator
% n: order of approximation (CFE expansion);
% r: the fractional order
% ts: sampling time in sec.
% Limit: max. order of CFE is 9
% this program is available at 
% http://www.crosswinds.net/~yqchen/software/foc/index.html

function [b,a]=tustin_cfe(n,r,ts)
constant0=(2 /ts)^r;
if (n>9) 
    n=9;
    disp('Max. term is clipped to 9');
end

if (n==1);
b(1)= 1/r;
b(2)= -1;
a(1)= 1/r;
a(2)= 1;
elseif (n==2) 
b(1)= 3/(r^2-1);
b(2)= -3*r/(r^2-1);
b(3)= 1;
a(1)= 3/(r^2-1);
a(2)= 3*r/(r^2-1);
a(3)= 1;

elseif (n==3)
b(1)=15/(r^2-4)/r;
b(2)=-15/(r^2-4);
b(3)=-(9-6*r^2)/(r^2-4)/r;
b(4)=-(-4*r+r^3)/(r^2-4)/r;
a(1)=15/(r^2-4)/r;
a(2)=15/(r^2-4);
a(3)=+(-9+6*r^2)/(r^2-4)/r;
a(4)=(-4*r+r^3)/(r^2-4)/r;

elseif (n==4)
b(1)=105/(r^2-1)/(r^2-9);
b(2)=-105*r/(r^2-1)/(r^2-9);
b(3)=+(45*r^2-90)/(r^2-1)/(r^2-9);
b(4)=+(-10*r^3+55*r)/(r^2-1)/(r^2-9);
b(5)=(r^4-10*r^2+9)/(r^2-1)/(r^2-9);
a(1)=105/(r^2-1)/(r^2-9);
a(2)=105*r/(r^2-1)/(r^2-9);
a(3)=+(45*r^2-90)/(r^2-1)/(r^2-9);
a(4)=+(10*r^3-55*r)/(r^2-1)/(r^2-9);
a(5)=(r^4-10*r^2+9)/(r^2-1)/(r^2-9);

elseif (n==5)
b(1)=945/(r^2-4)/(r^2-16)/r;
b(2)=-945/(r^2-4)/(r^2-16);
b(3)=-(-420*r^2+1050)/(r^2-4)/(r^2-16)/r;
b(4)=-(105*r^3-735*r)/(r^2-4)/(r^2-16)/r;
b(5)=-(195*r^2-225-15*r^4)/(r^2-4)/(r^2-16)/r;
b(6)=-(r^5+64*r-20*r^3)/(r^2-4)/(r^2-16)/r;
a(1)=945/(r^2-4)/(r^2-16)/r;
a(2)=945/(r^2-4)/(r^2-16);
a(3)=+(420*r^2-1050)/(r^2-4)/(r^2-16)/r;
a(4)=+(105*r^3-735*r)/(r^2-4)/(r^2-16)/r;
a(5)=+(-195*r^2+225+15*r^4)/(r^2-4)/(r^2-16)/r;
a(6)=(r^5+64*r-20*r^3)/(r^2-4)/(r^2-16)/r;

elseif (n==6)
b(1)=10395/(r^2-1)/(r^2-9)/(r^2-25);
b(2)=-10395*r/(r^2-1)/(r^2-9)/(r^2-25);
b(3)=+(-14175+4725*r^2)/(r^2-1)/(r^2-9)/(r^2-25);
b(4) = (10710*r-1260*r^3)/(r^2-1)/(r^2-9)/(r^2-25);
b(5) = (4725+210*r^4-3360*r^2)/(r^2-1)/(r^2-9)/(r^2-25);
b(6) = (525*r^3-2079*r-21*r^5)/(r^2-1)/(r^2-9)/(r^2-25);
b(7)=(r^6-35*r^4+259*r^2-225)/(r^2-1)/(r^2-9)/(r^2-25);
a(1)=10395/(r^2-1)/(r^2-9)/(r^2-25);
a(2)= 10395*r/(r^2-1)/(r^2-9)/(r^2-25);
a(3) = (-14175+4725*r^2)/(r^2-1)/(r^2-9)/(r^2-25);
a(4) = (-10710*r+1260*r^3)/(r^2-1)/(r^2-9)/(r^2-25);
a(5) = (4725+210*r^4-3360*r^2)/(r^2-1)/(r^2-9)/(r^2-25);
a(6) = (-525*r^3+2079*r+21*r^5)/(r^2-1)/(r^2-9)/(r^2-25);
a(7)=(r^6-35*r^4+259*r^2-225)/(r^2-1)/(r^2-9)/(r^2-25);

elseif (n==7)
b(1)=135135/(r^2-4)/(r^2-16)/(r^2-36)/r;
b(2)=-135135/(r^2-4)/(r^2-16)/(r^2-36);
b(3)=-(-62370*r^2+218295)/(r^2-4)/(r^2-16)/(r^2-36)/r;
b(4)=-(-173250*r+17325*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/r;
b(5)=-(-99225+59850*r^2-3150*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/r;
b(6)=-(53487*r+378*r^5-11340*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/r;
b(7)=-(1190*r^4+11025-28*r^6-10612*r^2)/(r^2-4)/(r^2-16)/(r^2-36)/r;
b(8)=-(-2304*r+r^7+784*r^3-56*r^5)/(r^2-4)/(r^2-16)/(r^2-36)/r;
a(1)=135135/(r^2-4)/(r^2-16)/(r^2-36)/r;
a(2)=135135/(r^2-4)/(r^2-16)/(r^2-36);
a(3) = (62370*r^2-218295)/(r^2-4)/(r^2-16)/(r^2-36)/r;
a(4) = (-173250*r+17325*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/r;
a(5) = (99225-59850*r^2+3150*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/r;
a(6) = (53487*r+378*r^5-11340*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/r;
a(7) = (-1190*r^4-11025+28*r^6+10612*r^2)/(r^2-4)/(r^2-16)/(r^2-36)/r;
a(8)=(-2304*r+r^7+784*r^3-56*r^5)/(r^2-4)/(r^2-16)/(r^2-36)/r;

elseif (n==8)
b(1)=2027025/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(2)=-2027025*r/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(3) = (945945*r^2-3783780)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(4) = (3108105*r-270270*r^3)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(5) = (2182950-1143450*r^2+51975*r^4)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(6) = (-1327095*r+242550*r^3-6930*r^5)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(7) = (-396900+328545*r^2+630*r^6-31500*r^4)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(8) = (-36*r^7+2394*r^5-39564*r^3+136431*r)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
b(9)=(r^8-84*r^6+1974*r^4-12916*r^2+11025)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(1)=2027025/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(2)=2027025*r/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(3) = (945945*r^2-3783780)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(4) = (-3108105*r+270270*r^3)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(5) = (2182950-1143450*r^2+51975*r^4)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(6) = (1327095*r-242550*r^3+6930*r^5)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(7) = (-396900+328545*r^2+630*r^6-31500*r^4)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(8) = (36*r^7-2394*r^5+39564*r^3-136431*r)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
a(9)=(r^8-84*r^6+1974*r^4-12916*r^2+11025)/(r^2-1)/(r^2-9)/(r^2-25)/(r^2-49);
elseif (n==9)
    b(1)=34459425/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(2)=-34459425/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64);
    b(3)=-(72972900-16216200*r^2)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(4)=-(-61486425*r+4729725*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(5)=-(-51081030+23648625*r^2-945945*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(6)=-(33648615*r+135135*r^5-5405400*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(7)=-(-13860*r^6-9514890*r^2+13097700+796950*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(8)=-(-5742495*r+1451835*r^3+990*r^7-76230*r^5)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(9)=-(-45*r^8+909765*r^2-893025+4410*r^6-120330*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    b(10)=-(4368*r^5-52480*r^3+r^9+147456*r-120*r^7)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(1)=34459425/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(2)=34459425/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64);
    a(3) = (-72972900+16216200*r^2)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(4) = (-61486425*r+4729725*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(5) = (51081030-23648625*r^2+945945*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(6) = (33648615*r+135135*r^5-5405400*r^3)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(7) = (13860*r^6+9514890*r^2-13097700-796950*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(8) = (-5742495*r+1451835*r^3+990*r^7-76230*r^5)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(9) = (45*r^8-909765*r^2+893025-4410*r^6+120330*r^4)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
    a(10)=(4368*r^5-52480*r^3+r^9+147456*r-120*r^7)/(r^2-4)/(r^2-16)/(r^2-36)/(r^2-64)/r;
%else
end
b=b.*constant0;
return