main.js (5830B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2019 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 * 18 * 19 * ## Notice 20 * 21 * The original Julia code and copyright notice are from the [Julia library]{@link https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/src/ellip.jl}. The implementation has been modified for JavaScript. 22 * 23 * ```text 24 * The MIT License (MIT) 25 * 26 * Copyright (c) 2017 Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and others: 27 * 28 * https://github.com/JuliaMath/SpecialFunctions.jl/graphs/contributors 29 * 30 * Permission is hereby granted, free of charge, to any person obtaining a copy 31 * of this software and associated documentation files (the "Software"), to deal 32 * in the Software without restriction, including without limitation the rights 33 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 34 * copies of the Software, and to permit persons to whom the Software is 35 * furnished to do so, subject to the following conditions: 36 * 37 * The above copyright notice and this permission notice shall be included in all 38 * copies or substantial portions of the Software. 39 * 40 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 41 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 42 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 43 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 44 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 45 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 46 * SOFTWARE. 47 * ``` 48 */ 49 50 'use strict'; 51 52 // MODULES // 53 54 var sqrt = require( './../../../../base/special/sqrt' ); 55 var ln = require( './../../../../base/special/ln' ); 56 var PINF = require( '@stdlib/constants/float64/pinf' ); 57 var HALF_PI = require( '@stdlib/constants/float64/half-pi' ); 58 var poly1 = require( './poly_p1.js' ); 59 var poly2 = require( './poly_p2.js' ); 60 var poly3 = require( './poly_p3.js' ); 61 var poly4 = require( './poly_p4.js' ); 62 var poly5 = require( './poly_p5.js' ); 63 var poly6 = require( './poly_p6.js' ); 64 var poly7 = require( './poly_p7.js' ); 65 var poly8 = require( './poly_p8.js' ); 66 var poly9 = require( './poly_p9.js' ); 67 var poly10 = require( './poly_p10.js' ); 68 var poly11 = require( './poly_p11.js' ); 69 var poly12 = require( './poly_p12.js' ); 70 71 72 // VARIABLES // 73 74 var ONE_DIV_PI = 0.3183098861837907; 75 76 77 // MAIN // 78 79 /** 80 * Computes the complete elliptic integral of the first kind. 81 * 82 * ## Method 83 * 84 * - The function computes the complete elliptic integral of the first kind in terms of parameter \\( m \\), instead of the elliptic modulus \\( k \\). 85 * 86 * ```tex 87 * K(m) = \int_0^{\pi/2} \frac{1}{\sqrt{1 - m sin^2\theta}} d\theta 88 * ``` 89 * 90 * - The function uses a piecewise approximation polynomial as given in Fukushima (2009). 91 * 92 * - For \\( m < 0 \\), the implementation follows Fukushima (2015). Namely, we use Equation 17.4.17 from the _Handbook of Mathematical Functions_ (Abramowitz and Stegun) to compute the function for \\( m < 0 \\) in terms of the piecewise polynomial representation of \\( m > 0 )). 93 * 94 * ```tex 95 * F(\phi|-m) = (1+m)^(-1/2) K(m/(1+m)) - (1+m)^(-1/2) F(\pi/2-\phi|m/(1+m)) 96 * ``` 97 * 98 * Since \\( K(m) \\) is equivalent to \\( F(\phi|m) \\), the above reduces to 99 * 100 * ```tex 101 * F(\phi|-m) = (1+m)^(-1/2) K(m/(1+m)) 102 * ``` 103 * 104 * ## References 105 * 106 * - Fukushima, Toshio. 2009. "Fast computation of complete elliptic integrals and Jacobian elliptic functions." _Celestial Mechanics and Dynamical Astronomy_ 105 (4): 305. doi:[10.1007/s10569-009-9228-z](https://doi.org/10.1007/s10569-009-9228-z). 107 * - Fukushima, Toshio. 2015. "Precise and fast computation of complete elliptic integrals by piecewise minimax rational function approximation." _Journal of Computational and Applied Mathematics_ 282 (July): 71–76. doi:[10.1016/j.cam.2014.12.038](https://doi.org/10.1016/j.cam.2014.12.038). 108 * 109 * @param {number} m - input value 110 * @returns {number} evaluated elliptic integral 111 * 112 * @example 113 * var v = ellipk( 0.5 ); 114 * // returns ~1.854 115 * 116 * v = ellipk( 2.0 ); 117 * // returns NaN 118 * 119 * v = ellipk( -1.0 ); 120 * // returns ~1.311 121 * 122 * v = ellipk( Infinity ); 123 * // returns NaN 124 * 125 * v = ellipk( -Infinity ); 126 * // returns NaN 127 * 128 * v = ellipk( NaN ); 129 * // returns NaN 130 */ 131 function ellipk( m ) { 132 var FLG; 133 var kdm; 134 var td; 135 var qd; 136 var t; 137 var x; 138 139 x = m; 140 if ( m < 0.0 ) { 141 x = m / ( m - 1.0 ); 142 FLG = true; 143 } 144 if ( x === 0.0 ) { 145 return HALF_PI; 146 } 147 if ( x === 1.0 ) { 148 return PINF; 149 } 150 if ( x > 1.0 ) { 151 return NaN; 152 } 153 if ( x < 0.1 ) { 154 t = poly1( x - 0.05 ); 155 } else if ( x < 0.2 ) { 156 t = poly2( x - 0.15 ); 157 } else if ( x < 0.3 ) { 158 t = poly3( x - 0.25 ); 159 } else if ( x < 0.4 ) { 160 t = poly4( x - 0.35 ); 161 } else if ( x < 0.5 ) { 162 t = poly5( x - 0.45 ); 163 } else if ( x < 0.6 ) { 164 t = poly6( x - 0.55 ); 165 } else if ( x < 0.7 ) { 166 t = poly7( x - 0.65 ); 167 } else if ( x < 0.8 ) { 168 t = poly8( x - 0.75 ); 169 } else if ( x < 0.85 ) { 170 t = poly9( x - 0.825 ); 171 } else if ( x < 0.9 ) { 172 t = poly10( x - 0.875 ); 173 } else { 174 td = 1.0 - x; 175 qd = poly11( td ); 176 kdm = poly12( td - 0.05 ); 177 t = -ln( qd ) * ( kdm * ONE_DIV_PI ); 178 } 179 if ( FLG ) { 180 // Complete the transformation mentioned above for m < 0: 181 return t / sqrt( 1.0 - m ); 182 } 183 return t; 184 } 185 186 187 // EXPORTS // 188 189 module.exports = ellipk;