188ef7ba56c77bc64de6cf371ca180c6418eb75d
[scilab.git] / scilab / modules / cacsd / macros / repfreq.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 1984 - 2011 - INRIA - Serge STEER
3 //
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function [frq,rep,splitf]=repfreq(sys,fmin,fmax,pas)
14
15     pas_def="auto";
16     l10=log(10);
17     [lhs,rhs]=argn(0)
18     //discretization
19     if and(typeof(sys)<>[ "rational" "state-space" "zpk"]) then
20         args=["sys","fmin","fmax","pas"];
21         ierr=execstr("[frq,rep,splitf]=%"+overloadname(sys)+"_repfreq("+strcat(args(1:rhs),",")+")","errcatch")
22         if ierr<>0 then
23             error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear dynamical system expected.\n"),"repfreq",1))
24         end
25         return
26     end
27     if typeof(sys)=="zpk" then sys=zpk2tf(sys);end
28     dom=sys.dt
29     if dom==[]|dom==0 then error(96,1),end
30     if dom=="d" then dom=1;end
31
32     select  rhs
33     case 1 then
34         pas=pas_def
35         if dom=="c" then fmax=1.d3; else fmax=1/(2*dom),end
36         fmin=0
37     case 2 then
38         if type(fmin)==10 then
39             rhs=1
40             pas=pas_def
41             if dom=="c" then fmax=1.d3; else fmax=1/(2*dom),end
42             if fmin=="sym" then
43                 fmin=-fmax
44             else
45                 fmin=0
46             end
47         else
48             frq=fmin
49         end
50     case 3 then
51         pas=pas_def
52     case 4 then ,
53     else
54         error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"), "repfreq",1,4))
55     end;
56     splitf=1
57     if rhs<>2 then
58         if fmin==[] then
59             fmin=0,
60         end
61         if fmax==[]|fmax==%inf then
62             if dom=="c" then
63                 fmax=1.d3;
64             else
65                 fmax=1/(2*dom);
66             end
67         end
68
69         if type(pas)==1 then
70             splitf=1
71             eps=1.e-14
72             if fmin<0&fmax>=0 then
73                 frq=- [exp(l10*((log(eps)/l10):pas:(log(-fmin)/l10))) -fmin];
74                 if fmax>eps then
75                     frq1=[exp(l10*((log(eps)/l10):pas:(log(fmax)/l10))) fmax];
76                     frq=[frq($:-1:1) frq1]
77                 else
78                     frq=frq($:-1:1);
79                 end
80             elseif fmin<0&fmax<0 then
81                 frq= [exp(l10*((log(-fmax)/l10):pas:(log(-fmin)/l10))) -fmin];
82                 frq=-frq($:-1:1);
83             elseif fmin >= fmax then
84                 error(msprintf(gettext("%s: Wrong value for input arguments #%d and #%d: %s < %s expected.\n"),..
85                 "repfreq",2,3,"fmin","fmax"));
86             else
87                 fmin=max(eps,fmin);
88                 frq=[exp(l10*((log(fmin)/l10):pas:(log(fmax)/l10))) fmax];
89             end
90         else
91             [frq,bnds,splitf]=calfrq(sys,fmin,fmax)
92         end;
93     end
94     //
95     typ=sys(1)
96     select typ(1)
97     case "r" then
98         [n,d]=sys(["num","den"]),
99         [mn,nn]=size(n)
100         if nn<>1 then error(95,1),end
101         if dom=="c" then
102             rep=freq(n,d,2*%pi*%i*frq),
103         else
104             rep=freq(n,d,exp(2*%pi*%i*dom*frq)),
105         end;
106     case "lss" then
107         [a,b,c,d]=abcd(sys)
108         [mn,nn]=size(d)
109         if nn<>1 then error(95,1),end
110
111         if dom=="c" then
112             if a==[] then
113                 rep=(d*2*%pi*%i)*frq
114             else
115                 rep=freq(a,b,c,d,2*%pi*%i*frq)
116             end
117         else
118             if a==[] then
119                 rep=d*exp(2*%pi*%i*dom*frq)
120             else
121                 rep=freq(a,b,c,d,exp(2*%pi*%i*dom*frq))
122             end
123         end;
124     else error(97,1),
125     end;
126     //representation
127     if lhs==1 then frq=rep,end
128 endfunction