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