function sysdfod=dfod1(n,T,a,r)
% sysdfod=dfod1(n,T,a,r): new digital fractional order differentiator
%                         and integrator 
%
% Output: =>
% Discrete system in the form of the rational approximation - IIR filter
% obtained by continued fraction expansion of a new generating function.
%
% Inputs: <=
% n: order of truncation n:=(1 2 3 4 5)
% T: sampling period in [sec]
% a: weighting factor between Euler and Tustin rules ( 0 <= a <= 1 )
%    0 - Euler rule, 1 - Tustin rule, 1/7 - Al-Alaoui rule, etc.
% r: approximated fractional order (s^r), r is generally real number
%
% Author: Dr. Ivo Petras (ivo.petras@mail.com)
% URL: http://ivopetras.tripod.com/
%
% Note: This approach is based on the original Al-Alaoui's work.
%
% Copyright (c), 2003.
%
k = ((1+a)/T)^r;
if (n>5) 
    n=5;
    disp('Maximal term is clipped to 5!');
end
if (n==1)
    b0=2/(a*r+a+r-1);
    b1=-(r+a*r-a+1)/(a*r+a+r-1);
    a0=2/(a*r+a+r-1);
    a1=1;
    sys1=tf(k*[b0 b1],[a0 a1],T,'Variable', 'z^-1');
    sysdfod=sys1;
elseif (n==2)
    b0=12/(a^2*r^2+3*a^2*r+2*a*r^2+2*a^2+r^2-8*a-3*r+2);
    b1=(12*a-6*r-6*a*r-12)/(a^2*r^2+3*a^2*r+2*a*r^2+2*a^2+r^2-8*a-3*r+2);
    b2=(a^2*r^2-3*a^2*r+2*a*r^2+r^2-8*a+3*r+2*a^2+2)/(a^2*r^2+3*a^2*r+2*a*r^2+2*a^2+r^2-8*a-3*r+2);
    a0=12/(a^2*r^2+3*a^2*r+2*a*r^2+2*a^2+r^2-8*a-3*r+2);
    a1=(12*a+6*r+6*a*r-12)/(a^2*r^2+3*a^2*r+2*a*r^2+2*a^2+r^2-8*a-3*r+2)
    a2=1;
    sys2=tf(k*[b0 b1 b2],[a0 a1 a2],T,'Variable', 'z^-1');
    sysdfod=sys2;
elseif (n==3)
    b0=120/(a^3*r^3+6*a^3*r^2+3*a^2*r^3+11*a^3*r+6*a^2*r^2+3*a*r^3+6*a^3-27*a^2*r-6*a*r^2+r^3-54*a^2-27*a*r-6*r^2+54*a+11*r-6);
    b1=-(60*r+180-180*a+60*a*r)/(a^3*r^3+6*a^3*r^2+3*a^2*r^3+11*a^3*r+6*a^2*r^2+3*a*r^3+6*a^3-27*a^2*r-6*a*r^2+r^3-54*a^2-27*a*r-6*r^2+54*a+11*r-6);
    b2=-(-72-72*a^2-24*a*r^2+216*a-60*r-12*r^2-12*a^2*r^2+60*a^2*r)/(a^3*r^3+6*a^3*r^2+3*a^2*r^3+11*a^3*r+6*a^2*r^2+3*a*r^3+6*a^3-27*a^2*r-6*a*r^2+r^3-54*a^2-27*a*r-6*r^2+54*a+11*r-6);
    b3=-(11*r+6*r^2+r^3+6*a*r^2-6*a^2*r^2-27*a*r-6*a^3*r^2+3*a^2*r^3+a^3*r^3+3*a*r^3+11*a^3*r-27*a^2*r-6*a^3-54*a+6+54*a^2)/(a^3*r^3+6*a^3*r^2+3*a^2*r^3+11*a^3*r+6*a^2*r^2+3*a*r^3+6*a^3-27*a^2*r-6*a*r^2+r^3-54*a^2-27*a*r-6*r^2+54*a+11*r-6);
    a0=120/(a^3*r^3+6*a^3*r^2+3*a^2*r^3+11*a^3*r+6*a^2*r^2+3*a*r^3+6*a^3-27*a^2*r-6*a*r^2+r^3-54*a^2-27*a*r-6*r^2+54*a+11*r-6);
    a1=(60*r-180+180*a+60*a*r)/(a^3*r^3+6*a^3*r^2+3*a^2*r^3+11*a^3*r+6*a^2*r^2+3*a*r^3+6*a^3-27*a^2*r-6*a*r^2+r^3-54*a^2-27*a*r-6*r^2+54*a+11*r-6);
    a2=(72+72*a^2+24*a*r^2-216*a-60*r+12*r^2+12*a^2*r^2+60*a^2*r)/(a^3*r^3+6*a^3*r^2+3*a^2*r^3+11*a^3*r+6*a^2*r^2+3*a*r^3+6*a^3-27*a^2*r-6*a*r^2+r^3-54*a^2-27*a*r-6*r^2+54*a+11*r-6);
    a3=1;
    sys3=tf(k*[b0 b1 b2 b3],[a0 a1 a2 a3],T,'Variable', 'z^-1');
    sysdfod=sys3;
elseif (n==4)
    b0=1680/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    b1=(3360*a-840*r-840*a*r-3360)/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    b2=(-5760*a+2160*a^2+2160+360*a*r^2+1260*r+180*r^2-1260*a^2*r+180*a^2*r^2)/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    b3=(-520*r+960*a^2*r-480-60*a^2*r^3+480*a^3-60*a*r^3+180*a^3*r^2-520*a^3*r+180*a^2*r^2+2880*a-180*r^2-180*a*r^2-20*r^3+960*a*r-20*a^3*r^3-2880*a^2)/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    b4=(10*r^3-40*a*r^2+r^4-150*a^2*r^2+35*r^2-384*a^3-384*a+50*r-320*a*r+24+20*a*r^3+a^4*r^4-10*a^4*r^3+4*a^3*r^4+35*a^4*r^2-20*a^3*r^3+6*a^2*r^4-50*a^4*r-40*a^3*r^2+4*a*r^4+864*a^2+24*a^4+320*a^3*r)/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    a0=1680/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    a1=(3360*a+840*r+840*a*r-3360)/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    a2=(-5760*a+2160*a^2+2160+360*a*r^2-1260*r+180*r^2+1260*a^2*r+180*a^2*r^2)/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    a3=(520*r-960*a^2*r-480+60*a^2*r^3+480*a^3+60*a*r^3+180*a^3*r^2+520*a^3*r+180*a^2*r^2+2880*a-180*r^2-180*a*r^2+20*r^3-960*a*r+20*a^3*r^3-2880*a^2)/(a^4*r^4+10*a^4*r^3+4*a^3*r^4+35*a^4*r^2+20*a^3*r^3+6*a^2*r^4+50*a^4*r-40*a^3*r^2+4*a*r^4+24*a^4-320*a^3*r-150*a^2*r^2-20*a*r^3+r^4-384*a^3-40*a*r^2-10*r^3+864*a^2+320*a*r+35*r^2-384*a-50*r+24);
    a4=1;
    sys4=tf(k*[b0 b1 b2 b3 b4],[a0 a1 a2 a3 a4],T,'Variable', 'z^-1');
    sysdfod=sys4;
elseif (n==5)
    b0=30240/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    b1=-(-75600*a+15120*r+15120*a*r+75600)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    b2=-(168000*a-30240*r-67200-3360*a^2*r^2-3360*r^2-67200*a^2-6720*a*r^2+30240*a^2*r)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    b3=-(5040*a*r^2-31500*a^2*r-5040*a^2*r^2-25200*a^3+19740*a^3*r+5040*r^2+1260*a*r^3+126000*a^2-31500*a*r+420*a^3*r^3+19740*r+420*r^3+1260*a^2*r^3+25200-5040*a^3*r^2-126000*a)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    b4=-(36000*a-120*a*r^4-120*a^3*r^4-72000*a^2+420*a^4*r^3+1560*a^3*r^2-3600-4620*r+1560*a*r^2-420*r^3-180*a^2*r^4-3600*a^4-21000*a^3*r+7380*a^2*r^2-840*a*r^3-30*r^4+21000*a*r-2130*r^2-2130*a^4*r^2+840*a^3*r^3-30*a^4*r^4+36000*a^3+4620*a^4*r)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    b5=-(120-3000*a+274*r+12000*a^2-12000*a^3+3000*a^4-120*a^5+4000*a^2*r+225*r^2-1005*a*r^2+5*a^4*r^3+1005*a^4*r^2-3250*a^4*r+1230*a^3*r^2+5*a*r^3-410*a^2*r^3-410*a^3*r^3-1230*a^2*r^2+4000*a^3*r+15*r^4+r^5-45*a^4*r^4+45*a*r^4+30*a^2*r^4-30*a^3*r^4+85*a^5*r^3-225*a^5*r^2+274*a^5*r-3250*a*r+85*r^3-15*a^5*r^4+5*a^4*r^5+5*a*r^5+10*a^2*r^5+10*a^3*r^5+a^5*r^5)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    a0=30240/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    a1=(75600*a+15120*r+15120*a*r-75600)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    a2=(-168000*a-30240*r+67200+3360*a^2*r^2+3360*r^2+67200*a^2+6720*a*r^2+30240*a^2*r)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    a3=(-5040*a*r^2-31500*a^2*r+5040*a^2*r^2+25200*a^3+19740*a^3*r-5040*r^2+1260*a*r^3-126000*a^2-31500*a*r+420*a^3*r^3+19740*r+420*r^3+1260*a^2*r^3-25200+5040*a^3*r^2+126000*a)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    a4=(-36000*a+120*a*r^4+120*a^3*r^4+72000*a^2+420*a^4*r^3-1560*a^3*r^2+3600-4620*r-1560*a*r^2-420*r^3+180*a^2*r^4+3600*a^4-21000*a^3*r-7380*a^2*r^2-840*a*r^3+30*r^4+21000*a*r+2130*r^2+2130*a^4*r^2+840*a^3*r^3+30*a^4*r^4-36000*a^3+4620*a^4*r)/(a^5*r^5+15*a^5*r^4+5*a^4*r^5+85*a^5*r^3+45*a^4*r^4+10*a^3*r^5+225*a^5*r^2+5*a^4*r^3+30*a^3*r^4+10*a^2*r^5+274*a^5*r-1005*a^4*r^2-410*a^3*r^3-30*a^2*r^4+5*a*r^5+120*a^5-3250*a^4*r-1230*a^3*r^2-410*a^2*r^3-45*a*r^4+r^5-3000*a^4+4000*a^3*r+1230*a^2*r^2+5*a*r^3-15*r^4+12000*a^3+4000*a^2*r+1005*a*r^2+85*r^3-12000*a^2-3250*a*r-225*r^2+3000*a+274*r-120);
    a5=1;
    sys5=tf(k*[b0 b1 b2 b3 b4 b5],[a0 a1 a2 a3 a4 a5],T,'Variable', 'z^-1');
    sysdfod=sys5;
end          
%