1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2//
3// Copyright (C) 2019 - Samuel GOUGEON - Le Mans Université
4//
5// This file is hereby licensed under the terms of the GNU GPL v2.0,
6// pursuant to article 5.3.4 of the CeCILL v.2.1.
7// This file was originally licensed under the terms of the CeCILL v2.1,
8// and continues to be available under such terms.
9// For more information, see the COPYING file which you should have received
10// along with this program.
11
12function varargout = ellipj(varargin)
13    //           sn = ellipj(x, m)  // Replaces %sn(x,m)
14    //     [sn, cn] = ellipj(x, m)  // Returns the cn function in addition to sn
15    // [sn, cn, dn] = ellipj(x, m)  // Returns dn function in addition to sn and cn
16    //           st = ellipj(funNames, x, m)
17    //               // Returns the Jacobi functions named in the funNames
18    //               //  vector of strings among "sn" "cn" "dn" "ns" "nc" "nd".
19    //               //  Values are returned through the st structure whose fields
20    //               //  are the elements of funNames.
21    //               //  Later, this syntax could be easily extendable to any
22    //               //  of the other non trivial "cd" "cs" "dc" "ds" "sc"
23    //               //  and "sd" Jacobi functions.
24
25    fname = "ellipj"
26    lhs = argn(1)
27    funs = []   // names of Jacobi functions to be returned
28    shift = 0   // input argument shift, for error messages about x or m
29
30    // CHECKING ARGUMENTS
31    // ==================
32    in = varargin
33    sizin = length(varargin)
34    if sizin < 2 | sizin > 3 then
35        msg = _("%s: Wrong number of input arguments: %d or %d expected.\n")
36        error(msprintf(msg, fname, 2, 3))
37    end
38    if sizin==3 then
39        [funs, k]= unique(in(1))
40        funs = in(1)(gsort(k,"g","i"))
41        in(1) = null()
42        if lhs > 1
43            msg = _("%s: Wrong number of output arguments: %d expected.\n")
44            error(msprintf(msg, fname, 1))
45        end
46        shift = 1
47        if type(funs) <> 10
48            msg = _("%s: Argument #%d: Text expected.\n")
49            error(msprintf(msg, fname, 1))
50        end
51        if length(grep(["cn" "dn" "sn" "nc" "nd" "ns"], funs)) <> size(funs,"*")
52            msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
53            error(msprintf(msg, fname, 1, "''cn'' ''dn'' ''sn'' ''nc'' ''nd'' ''ns''"))
54        end
55    else
56        if lhs > 3
57            msg = _("%s: Wrong number of output arguments: %d to %d expected.\n")
58            error(msprintf(msg, fname, 1, 3))
59        end
60    end
61    // x
62    x = in(1)
63    if type(x) <> 1 then
64        msg = _("%s: Argument #%d: Decimal or complex number expected.\n")
65        error(msprintf(msg, fname, 1+shift))
66    end
67    sx = size(x)
68    x = x(:)
69    // m
70    m = in(2)
71    if type(m) <> 1 | ~isreal(m,0) then
72        msg = _("%s: Argument #%d: Decimal number expected.\n")
73        error(msprintf(msg, fname, 2+shift))
74    end
75    if or(m < 0) | or(m > 1) then
76        msg = _("%s: Argument #%d: Must be in the interval [%d, %d].\n")
77        error(msprintf(msg, fname, 2+shift, 0, 1))
78    end
79
80    // PROCESSING
81    // ==========
82    [n1,n2] = size(x);
83    n = n1*n2;
84    a = amell(real(x), sqrt(m));
85    s = sin(a);
86    c = cos(a);
87    d = sqrt(ones(n1,n2) - m*s.*s);
88    xReal = isreal(x,0)
89    if ~xReal then
90        m1 = 1 - m;
91        a1 = amell(imag(x), sqrt(m1));
92        s1 = sin(a1);
93        c1 = cos(a1);
94        d1 = sqrt(ones(n1,n2) - m1*s1.*s1);
95        N = (c1.*c1 + m*s.*s.*s1.*s1);
96    end
97    kinf = isinf(x)
98
99    // sn
100    // --
101    if funs==[] | grep(funs,["sn" "ns" "dn" "nd"]) <> [] then
102        if xReal
103            sn = s
104        else
105            sn = (s.*d1 + %i*c.*d.*s1.*c1) ./ N;
106        end
107        if ~isreal(sn) & isreal(sn,0) then
108            sn = real(sn)
109        end
110        sn(kinf) = %nan
111        sn = matrix(sn, sx)
112    end
113    // cn
114    // --
115    if lhs > 1 | grep(funs,["cn" "nc"]) <> [] then
116        if xReal
117            cn = c
118        else
119            cn = (c.*c1 - %i*s.*d.*s1.*d1) ./ N;
120        end
121        if ~isreal(cn) & isreal(cn,0) then
122            cn = real(cn)
123        end
124        cn(kinf) = %nan
125        cn = matrix(cn, sx)
126    end
127    // dn
128    // --
129    if lhs > 2 | grep(funs,["dn" "nd"]) <> [] then
130        dn = sqrt(1 - m*sn.*sn);
131        if ~isreal(dn) & isreal(dn,0) then
132            dn = real(dn)
133        end
134    end
135    // nc, nd, ns
136    // ----------
137    if funs <> [] then
138        if grep(funs, "nc") <> [], nc = 1 ./ cn, end
139        if grep(funs, "nd") <> [], nd = 1 ./ dn, end
140        if grep(funs, "ns") <> [], ns = 1 ./ sn, end
141    end
142
143    // OUTPUT
144    // ======
145    if funs <> [] then
146        st = struct()
147        for f = funs(:)'
148            execstr("st(f) = "+f)
149        end
150        varargout = list(st)
151    else
152        varargout = list()
153        varargout(1) = sn
154        if lhs > 1, varargout(2) = cn, end
155        if lhs > 2, varargout(3) = dn, end
156    end
157endfunction
158