main.js (5354B)
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 HALF_PI = require( '@stdlib/constants/float64/half-pi' ); 56 var ellipk = require( './../../../../base/special/ellipk' ); 57 var poly1 = require( './poly_p1.js' ); 58 var poly2 = require( './poly_p2.js' ); 59 var poly3 = require( './poly_p3.js' ); 60 var poly4 = require( './poly_p4.js' ); 61 var poly5 = require( './poly_p5.js' ); 62 var poly6 = require( './poly_p6.js' ); 63 var poly7 = require( './poly_p7.js' ); 64 var poly8 = require( './poly_p8.js' ); 65 var poly9 = require( './poly_p9.js' ); 66 var poly10 = require( './poly_p10.js' ); 67 var poly11 = require( './poly_p11.js' ); 68 var poly12 = require( './poly_p12.js' ); 69 70 71 // MAIN // 72 73 /** 74 * Computes the complete elliptic integral of the second kind. 75 * 76 * ## Method 77 * 78 * - The function computes the complete elliptic integral of the second kind in terms of parameter \\( m \\), instead of the elliptic modulus \\( k \\). 79 * 80 * ```tex 81 * E(m) = \int_0^{\pi/2} \sqrt{1 - m (\sin\theta)^2} d\theta 82 * ``` 83 * 84 * - The function uses a piecewise approximation polynomial as given in Fukushima (2009). 85 * 86 * - For \\( m < 0 \\), the implementation follows Fukushima (2015). 87 * 88 * ## References 89 * 90 * - 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). 91 * - 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). 92 * 93 * @param {number} m - input value 94 * @returns {number} evaluated elliptic integral 95 * 96 * @example 97 * var v = ellipe( 0.5 ); 98 * // returns ~1.351 99 * 100 * v = ellipe( -1.0 ); 101 * // returns ~1.910 102 * 103 * v = ellipe( 2.0 ); 104 * // returns NaN 105 * 106 * v = ellipe( Infinity ); 107 * // returns NaN 108 * 109 * v = ellipe( -Infinity ); 110 * // returns NaN 111 * 112 * v = ellipe( NaN ); 113 * // returns NaN 114 */ 115 function ellipe( m ) { 116 var FLG; 117 var kdm; 118 var edm; 119 var td; 120 var km; 121 var t; 122 var x; 123 124 x = m; 125 if ( m < 0.0 ) { 126 x = m / ( m - 1.0 ); 127 FLG = true; 128 } 129 if ( x === 0.0 ) { 130 return HALF_PI; 131 } 132 if ( x === 1.0 ) { 133 return 1.0; 134 } 135 if ( x > 1.0 ) { 136 return NaN; 137 } 138 if ( x < 0.1 ) { 139 t = poly1( x - 0.05 ); 140 } else if ( x < 0.2 ) { 141 t = poly2( x - 0.15 ); 142 } else if ( x < 0.3 ) { 143 t = poly3( x - 0.25 ); 144 } else if ( x < 0.4 ) { 145 t = poly4( x - 0.35 ); 146 } else if ( x < 0.5 ) { 147 t = poly5( x - 0.45 ); 148 } else if ( x < 0.6 ) { 149 t = poly6( x - 0.55 ); 150 } else if ( x < 0.7 ) { 151 t = poly7( x - 0.65 ); 152 } else if ( x < 0.8 ) { 153 t = poly8( x - 0.75 ); 154 } else if ( x < 0.85 ) { 155 t = poly9( x - 0.825 ); 156 } else if ( x < 0.9 ) { 157 t = poly10( x - 0.875 ); 158 } else { 159 td = 0.95 - x; 160 kdm = poly11(td); 161 edm = poly12(td); 162 km = ellipk( x ); 163 164 // To avoid precision loss near 1, we use Eq. 33 from Fukushima (2009): 165 t = ( HALF_PI + ( km * (kdm - edm) ) ) / kdm; 166 } 167 if ( FLG ) { 168 // Complete the transformation mentioned above for m < 0: 169 return t * sqrt( 1.0 - m ); 170 } 171 return t; 172 } 173 174 175 // EXPORTS // 176 177 module.exports = ellipe;