Press n or j to go to the next uncovered block, b, p or k for the previous block.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 | 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 7x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 24031x 6x 6x 6x 6x 6x 6x 24025x 24025x 24025x 24025x 24025x 24031x 12006x 12006x 24025x 24031x 1x 1x 1x 1x 1x 1x 24024x 24031x 7842x 7842x 7842x 127920x 127920x 127920x 7842x 7842x 127920x 7842x 7842x 7842x 7842x 125376x 125376x 125376x 7842x 7842x 125376x 7842x 7842x 7842x 7842x 7842x 7842x 7842x 124446x 124446x 124446x 7842x 7842x 124446x 7842x 7842x 7842x 7842x 129252x 129252x 129252x 7842x 7842x 129252x 7842x 7842x 7842x 24031x 16182x 16182x 666x 666x 16182x 7806x 7806x 16182x 16182x 9006x 16182x 7176x 7176x 7176x 7176x 12x 12x 7176x 12x 12x 7176x 12x 12x 7176x 7176x 7176x 12x 12x 7176x 12x 12x 7176x 7176x 16182x 16182x 16182x 16182x 16182x 16182x 16182x 278754x 278754x 278754x 278754x 278754x 16182x 16182x 9006x 9006x 9006x 9006x 110112x 110112x 110112x 110112x 9006x 9006x 9006x 9006x 110112x 110112x 110112x 110112x 9006x 9006x 9006x 9006x 9006x 16182x 7176x 7176x 7176x 7176x 7176x 7176x 7176x 80736x 80736x 80736x 80736x 7176x 7176x 7176x 7176x 80724x 80724x 80724x 80724x 7176x 7176x 7176x 7176x 7176x 7176x 16182x 24024x 24031x 7x 7x | /** * @license Apache-2.0 * * Copyright (c) 2025 The Stdlib Authors. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * * * ## Notice * * The original C code and copyright notice are from the [xsf library]{@link https://github.com/scipy/xsf/blob/main/include/xsf/airy.h}. The implementation has been modified for JavaScript. * * ```text * Copyright 2024, SciPy * * Use, modification and distribution are subject to the * BSD 3-Clause License (See accompanying file LICENSE or copy at * https://github.com/scipy/xsf/blob/main/LICENSE) * ``` */ /* eslint-disable max-statements */ 'use strict'; // MODULES // var SQRT3 = require( '@stdlib/constants/float64/sqrt-three' ); var floor = require( '@stdlib/math/base/special/floor' ); var PINF = require( '@stdlib/constants/float64/pinf' ); var sqrt = require( '@stdlib/math/base/special/sqrt' ); var sin = require( '@stdlib/math/base/special/sin' ); var cos = require( '@stdlib/math/base/special/cos' ); var exp = require( '@stdlib/math/base/special/exp' ); var abs = require( '@stdlib/math/base/special/abs' ); var PI = require( '@stdlib/constants/float64/pi' ); // VARIABLES // var C1 = 0.355028053887817; var C2 = 0.258819403792807; var EPS = 1.0e-15; var ck = new Array( 52 ); // eslint-disable-line stdlib/no-new-array var dk = new Array( 52 ); // eslint-disable-line stdlib/no-new-array /** * Evaluates the Airy functions of the first and second kind, Ai(x) and Bi(x), and their first derivatives, Ai'(x) and Bi'(x), and assigns to a provided output array. * * @param {number} x - input value * @param {Array} out - output array * @param {integer} stride - output array stride * @param {NonNegativeInteger} offset - output array index offset * @returns {Array} array containing Ai(x), Ai'(x), Bi(x), and Bi'(x) */ function airy(x, out, stride, offset) { // @eslint-disable-line max-statements var kmax; var km2; var sai; var sad; var sbi; var sbd; var xp1; var xr2; var xar; var xr1; var xcs; var xss; var ssa; var sda; var ssb; var sdb; var xa; var xq; var xm; var xf; var rp; var xe; var gx; var df; var dg; var fx; var km; var r; var k; if ( x === PINF ) { out[offset] = 0.0; out[offset + stride] = 0.0; out[offset + (2 * stride)] = PINF; out[offset + (3 * stride)] = PINF; return out; } km2 = 0; xa = abs(x); xq = sqrt(xa); xm = 8.0; if (x > 0.0) { xm = 5.0; } if (x === 0.0) { out[offset] = C1; out[offset + stride] = -C2; out[offset + (2 * stride)] = SQRT3 * C1; out[offset + (3 * stride)] = SQRT3 * C2; return out; } if (xa <= xm) { fx = 1.0; r = 1.0; for (k = 1; k <= 40; k++) { r = r * x / (3.0 * k) * x / ((3.0 * k) - 1.0) * x; fx += r; if (abs(r) < abs(fx) * EPS) { break; } } gx = x; r = x; for (k = 1; k <= 40; k++) { r = r * x / (3.0 * k) * x / ((3.0 * k) + 1.0) * x; gx += r; if (abs(r) < abs(gx) * EPS) { break; } } out[offset] = (C1 * fx) - (C2 * gx); out[offset + (2 * stride)] = SQRT3 * ((C1 * fx) + (C2 * gx)); df = 0.5 * x * x; r = df; for (k = 1; k <= 40; k++) { r = r * x / (3.0 * k) * x / ((3.0 * k) + 2.0) * x; df += r; if (abs(r) < abs(df) * EPS) { break; } } dg = 1.0; r = 1.0; for (k = 1; k <= 40; k++) { r = r * x / (3.0 * k) * x / ((3.0 * k) - 2.0) * x; dg += r; if (abs(r) < abs(dg) * EPS) { break; } } out[offset + stride] = (C1 * df) - (C2 * dg); out[offset + (3 * stride)] = SQRT3 * ((C1 * df) + (C2 * dg)); } else { km = floor(24.5 - xa); if (xa < 6.0) { km = 14; } if (xa > 15.0) { km = 10; } if (x > 0.0) { kmax = km; } else { /* Choose cutoffs so that the remainder term in asymptotic expansion is epsilon size. The X<0 branch needs to be fast in order to make AIRYZO efficient */ if (xa > 70.0) { km = 3; } if (xa > 500.0) { km = 2; } if (xa > 1000.0) { km = 1; } km2 = km; if (xa > 150.0) { km2 = 1; } if (xa > 3000.0) { km2 = 0; } kmax = (2 * km) + 1; } xe = xa * xq / 1.5; xr1 = 1.0 / xe; xar = 1.0 / xq; xf = sqrt(xar); rp = 0.5641895835477563; r = 1.0; for (k = 1; k <= kmax; k++) { r = r * ((6.0 * k) - 1.0) / 216.0 * ((6.0 * k) - 3.0) / k; r *= ((6.0 * k) - 5.0) / ((2.0 * k) - 1.0); ck[k - 1] = r; dk[k - 1] = -((6.0 * k) + 1.0) / ((6.0 * k) - 1.0) * r; } if (x > 0.0) { sai = 1.0; sad = 1.0; r = 1.0; for (k = 1; k <= km; k++) { r *= -xr1; sai += ck[k - 1] * r; sad += dk[k - 1] * r; } sbi = 1.0; sbd = 1.0; r = 1.0; for (k = 1; k <= km; k++) { r *= xr1; sbi += ck[k - 1] * r; sbd += dk[k - 1] * r; } xp1 = exp(-xe); out[offset] = 0.5 * rp * xf * xp1 * sai; out[offset + stride] = -0.5 * rp / xf * xp1 * sad; out[offset + (2*stride)] = rp * xf / xp1 * sbi; out[offset + (3*stride)] = rp / xf / xp1 * sbd; } else { xcs = cos(xe + (PI / 4.0)); xss = sin(xe + (PI / 4.0)); ssa = 1.0; sda = 1.0; r = 1.0; xr2 = 1.0 / (xe * xe); for (k = 1; k <= km; k++) { r *= -xr2; ssa += ck[(2 * k) - 1] * r; sda += dk[(2 * k) - 1] * r; } ssb = ck[0] * xr1; sdb = dk[0] * xr1; r = xr1; for (k = 1; k <= km2; k++) { r *= -xr2; ssb += ck[2 * k] * r; sdb += dk[2 * k] * r; } out[offset] = rp * xf * ((xss * ssa) - (xcs * ssb)); out[offset + stride] = -rp / xf * ((xcs * sda) + (xss * sdb)); out[offset + (2*stride)] = rp * xf * ((xcs * ssa) + (xss * ssb)); out[offset + (3*stride)] = rp / xf * ((xss * sda) - (xcs * sdb)); } } return out; } module.exports = airy; |