macrorugo.ts 21.3 KB
Newer Older
Mathias Chouet's avatar
Mathias Chouet committed
1
import { CalculatorType } from "../compute-node";
Dorchies David's avatar
Dorchies David committed
2
import { Dichotomie } from "../dichotomie";
3
import { MacrorugoCompound } from "../index";
4
import { ParamCalculability, ParamFamily } from "../param/param-definition";
5
import { ParamValueMode } from "../param/param-value-mode";
6
import { SessionSettings } from "../session_settings";
7
import { Message, MessageCode } from "../util/message";
Dorchies David's avatar
Dorchies David committed
8
9
import { Result } from "../util/result";
import { MacrorugoParams } from "./macrorugo_params";
10
import { MRCInclination } from "./mrc-inclination";
11
import { FishPass } from "../fish_pass";
Dorchies David's avatar
Dorchies David committed
12

13
14
15
export enum MacroRugoFlowType {
    EMERGENT,
    QUASI_EMERGENT,
16
    SUBMERGED
17
18
}

19
export class MacroRugo extends FishPass {
Dorchies David's avatar
Dorchies David committed
20

21
22
23
24
25
26
27
    private static readonly g = 9.81;

    /** nu: water kinematic viscosity  */
    private static readonly nu = 1E-6;
    // Water at 20 °C has a kinematic viscosity of about 10−6 m2·s−1
    // (https://en.wikipedia.org/wiki/Viscosity#Kinematic_viscosity,_%CE%BD)

28
    /** Ratio between the width (perpendicular to flow) and the length (parallel to flow) of a cell (-) */
29
    private static readonly fracAxAy = 1;
30
31

    /** Limit between emergent and submerged flow */
32
    private static readonly limitSubmerg = 1.1;
33

34
    /** Flag for submerged Flow */
35
    private flowType: MacroRugoFlowType;
36
37

    /** Velocity at the bed (m.s-1) */
Dorchies David's avatar
Dorchies David committed
38
39
    private u0: number;

40
    private _cache: { [key: string]: number };
Dorchies David's avatar
Dorchies David committed
41
42
43

    constructor(prms: MacrorugoParams, dbg: boolean = false) {
        super(prms, dbg);
44
        this._cache = {};
Mathias Chouet's avatar
Mathias Chouet committed
45
        this._calcType = CalculatorType.MacroRugo;
46
47
        this._defaultCalculatedParam = this.prms.Q;
        this.resetDefaultCalculatedParam();
Dorchies David's avatar
Dorchies David committed
48
49
50
51
52
53
54
55
56
    }

    /**
     * paramètres castés au bon type
     */
    get prms(): MacrorugoParams {
        return this._prms as MacrorugoParams;
    }

57
58
    /**
     * Calcul du débit total, de la cote amont ou aval ou d'un paramètre d'une structure
59
     * @param sVarCalc Nom du paramètre à calculer
60
61
     * @param rInit Valeur initiale
     */
62
    public Calc(sVarCalc: string, rInit?: number): Result {
Mathias Chouet's avatar
Mathias Chouet committed
63
64
65
66
        // Si on ne force pas à CALCUL, les tests à la 5e décimale ne passent plus (macrorugo.spec.ts)
        const originalValueMode = this.getParameter(sVarCalc).valueMode;
        this.getParameter(sVarCalc).valueMode = ParamValueMode.CALCUL;

67
        const r: Result = super.Calc(sVarCalc, rInit);
68
69
        this.getParameter(sVarCalc).valueMode = originalValueMode;

70
        // La largeur de la rampe est-elle adéquate par rapport à la largeur de motif ax ?
71
72
73
74
75
        // (si le parent est une MacroRugoCompound avec radier incliné, on ignore ces problèmes)
        if (
            this.parent === undefined
            || (
                this.parent instanceof MacrorugoCompound
76
                && this.parent.properties.getPropValue("inclinedApron") === MRCInclination.NOT_INCLINED
77
78
79
            )
        ) {
            const ax: number = this.prms.PBD.v / Math.sqrt(this.prms.C.v);
80
            if (this.prms.B.v < ax - 0.001) { // B < ax, with a little tolerance
81
82
83
                const m = new Message(MessageCode.WARNING_RAMP_WIDTH_LOWER_THAN_PATTERN_WIDTH);
                m.extraVar.pattern = ax;
                r.resultElement.log.add(m);
84
            } else if (this.prms.B.v % (ax / 2) > 0.001 && this.prms.B.v % (ax / 2) < ax / 2 - 0.001) {
85
86
87
88
89
90
                const m = new Message(MessageCode.WARNING_RAMP_WIDTH_NOT_MULTIPLE_OF_HALF_PATTERN_WIDTH);
                m.extraVar.halfPattern = ax / 2;
                m.extraVar.lower = Math.floor(this.prms.B.v / ax * 2) * (ax / 2);
                m.extraVar.higher = Math.ceil(this.prms.B.v / ax * 2) * (ax / 2);
                r.resultElement.log.add(m);
            }
91
92
        }

93
94
95
96
97
98
99
100
101
        // La concentration est-elle dans les valeurs admissibles 8-20% (#284)
        if (this.parent === undefined) {
            if(this.prms.C.V < 0.08 || this.prms.C.V > 0.2) {
                r.resultElement.log.add(
                    new Message(MessageCode.WARNING_MACRORUGO_CONCENTRATION_OUT_OF_BOUNDS)
                );
            }
        }

102
103
        // Ajout des résultats complémentaires
        // Cote de fond aval
104
        r.resultElement.values.ZF2 = this.prms.ZF1.v - this.prms.If.v * this.prms.L.v;
105
        // Vitesse débitante
106
107
108
109
110
        let resVdeb = this.prms.Q.V / this.prms.B.v / this.prms.Y.v;
        if (isNaN(resVdeb)) {
            resVdeb = 0;
        }
        r.resultElement.values.Vdeb = resVdeb;
111
112
113
114
115
116
117
118
119
120
        if (this.flowType !== MacroRugoFlowType.SUBMERGED) {
            // Froude
            r.resultElement.values.Vg =  r.resultElement.values.Vdeb / (1 - Math.sqrt(MacroRugo.fracAxAy * this.prms.C.v));
            let resFr = r.resultElement.values.Vg / Math.sqrt(MacroRugo.g * this.prms.Y.v);
            if (isNaN(resFr)) { // if Y == 0
                resFr = 0;
            }
            r.resultElement.values.Fr = resFr;
            // Vitesse maximale
            r.resultElement.values.Vmax = this.r * r.resultElement.values.Vg * this.CalcfFr(resVdeb);
121
        }
122
        // Puissance dissipée
123
        r.resultElement.values.PV = 1000 * MacroRugo.g * r.resultElement.values.Vdeb * this.prms.If.v;
124
125
        // Type d'écoulement
        if (this.prms.Y.v / this.prms.PBH.v < 1) {
126
            r.resultElement.values.ENUM_MacroRugoFlowType = MacroRugoFlowType.EMERGENT;
127
        } else if (this.prms.Y.v / this.prms.PBH.v < MacroRugo.limitSubmerg) {
128
            r.resultElement.values.ENUM_MacroRugoFlowType = MacroRugoFlowType.QUASI_EMERGENT;
129
        } else {
130
            r.resultElement.values.ENUM_MacroRugoFlowType = MacroRugoFlowType.SUBMERGED;
131
132
        }
        // Vitesse et débit du guide technique
133
        /* let cQ: [number, number, number, number];
134
        let cV: [number, number, number];
135
        let hdk: number;
136
        if (r.resultElement.values.ENUM_MacroRugoFlowType === MacroRugoFlowType.SUBMERGED) {
137
            cQ = [0.955, 2.282, 0.466, -0.23];
138
            hdk = this.prms.PBH.v;
139
        } else {
140
            hdk = this.prms.PBD.v;
141
142
143
144
145
            if (Math.abs(this.prms.Cd0.v - 2) < 0.05) {
                cQ = [0.648, 1.084, 0.56, -0.456];
                cV = [3.35, 0.27, 0.53];
            } else {
                cQ = [0.815, 1.45, 0.557, -0.456];
Dorchies David's avatar
Dorchies David committed
146
                cV = [4.54, 0.32, 0.56];
147
            }
148
        } */
Dorchies David's avatar
Dorchies David committed
149
150
151
152
153
154
        if (this.prms.Y.v > 0 && this.prms.If.v > 0) {
            r.resultElement.values.Strickler =
                this.prms.Q.V / (Math.pow(this.prms.Y.v, 5 / 3) * this.prms.B.v * Math.pow(this.prms.If.v, 0.5));
        } else {
            r.resultElement.values.Strickler = 0;
        }
155
156
157
        return r;
    }

Dorchies David's avatar
Dorchies David committed
158
    public Equation(sVarCalc: string): Result {
159
        const Q = this.prms.Q.v;
Dorchies David's avatar
Dorchies David committed
160
        const q0 = Math.sqrt(2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v * (1 - (this.sigma * this.prms.C.v)) /
161
            (this.Cx * this.prms.C.v)) * this.prms.Y.v * this.prms.B.v;
162
163
        let r: Result;
        if (q0 > 0) {
164
            this.setFlowType();
165
            const dicho = new Dichotomie(this, "Q", false, this.resolveQ);
166
            r = dicho.Dichotomie(0, SessionSettings.precision, q0);
167
168
169
        } else {
            r = new Result(0);
        }
170
        this.prms.Q.v = Q;
Dorchies David's avatar
Dorchies David committed
171
        return r;
Dorchies David's avatar
Dorchies David committed
172
173
    }

174
175
176
177
178
179
180
181
182
183
184
185
    /**
     * paramétrage de la calculabilité des paramètres
     */
    protected setParametersCalculability() {
        this.prms.ZF1.calculability = ParamCalculability.FREE;
        this.prms.L.calculability = ParamCalculability.FREE;
        this.prms.Ks.calculability = ParamCalculability.FREE;
        this.prms.B.calculability = ParamCalculability.DICHO;
        this.prms.If.calculability = ParamCalculability.DICHO;
        this.prms.Q.calculability = ParamCalculability.EQUATION;
        this.prms.Y.calculability = ParamCalculability.DICHO;
        this.prms.C.calculability = ParamCalculability.DICHO;
186
        this.prms.PBD.calculability = ParamCalculability.FREE;
187
        this.prms.PBH.calculability = ParamCalculability.FREE;
188
        this.prms.Cd0.calculability = ParamCalculability.FREE;
189
190
    }

191
192
193
194
195
    protected setResultsUnits() {
        this._resultsUnits = {
            PV: "W/m³",
            Vdeb: "m/s",
            Vmax: "m/s",
196
            Vg: "m/s",
197
198
199
200
201
            ZF2: "m",
            Strickler: "SI"
        }
    }

Mathias Chouet's avatar
Mathias Chouet committed
202
    protected exposeResults() {
203
        this._resultsFamilies = {
204
            ZF2: ParamFamily.ELEVATIONS,
205
            Vdeb: ParamFamily.SPEEDS,
Mathias Chouet's avatar
Mathias Chouet committed
206
            Vmax: ParamFamily.SPEEDS,
207
            Vg: ParamFamily.SPEEDS,
Mathias Chouet's avatar
Mathias Chouet committed
208
209
            Fr: undefined,
            PV: undefined,
210
            Strickler: ParamFamily.STRICKLERS
211
212
213
        };
    }

214
215
216
217
218
219
220
221
222
223
224
    private setFlowType() {
        const hstar: number = this.prms.Y.v / this.prms.PBH.v;
        if (hstar > MacroRugo.limitSubmerg) {
            this.flowType = MacroRugoFlowType.SUBMERGED;
        } else if (hstar < 1) {
            this.flowType = MacroRugoFlowType.EMERGENT;
        } else {
            this.flowType = MacroRugoFlowType.QUASI_EMERGENT;
        }
    }

Dorchies David's avatar
Dorchies David committed
225
    /**
226
227
     * Equation from Cassan, L., Laurens, P., 2016. Design of emergent and submerged rock-ramp fish passes.
     * Knowledge & Management of Aquatic Ecosystems 45.
Dorchies David's avatar
Dorchies David committed
228
229
     * @param sVarCalc Variable à calculer
     */
Dorchies David's avatar
Dorchies David committed
230
    private resolveQ(): number {
Dorchies David's avatar
Dorchies David committed
231
232
        // Reset cached variables depending on Q (or not...)
        this._cache = {};
233
234

        /** adimensional water depth */
235
        const hstar: number = this.prms.Y.v / this.prms.PBH.v;
236

237
238
239
240
241
242
243
244
        switch (this.flowType) {
            case MacroRugoFlowType.SUBMERGED:
                return this.resolveQSubmerged();
            case MacroRugoFlowType.EMERGENT:
                return this.resolveQEmergent();
            case MacroRugoFlowType.QUASI_EMERGENT:
                const a = (hstar - 1) / (MacroRugo.limitSubmerg - 1);
                return (1 - a) * this.resolveQEmergent() + a * this.resolveQSubmerged();
Dorchies David's avatar
Dorchies David committed
245
        }
246
247
248
249
250
251
252
    }

    /**
     * Averaged velocity (m.s-1)
     */
    private get U0(): number {
        if (this._cache.U0 === undefined) {
Dorchies David's avatar
Dorchies David committed
253
            this._cache.U0 = this.prms.Q.v / this.prms.B.v / this.prms.Y.v;
254
255
        }
        return this._cache.U0;
Dorchies David's avatar
Dorchies David committed
256
257
258
    }

    private get CdChD(): number {
259
260
        if (this._cache.CdChD === undefined) {
            this._cache.CdChD = this.Cd * this.prms.C.v * this.prms.PBH.v / this.prms.PBD.v;
Dorchies David's avatar
Dorchies David committed
261
        }
262
        return this._cache.CdChD;
Dorchies David's avatar
Dorchies David committed
263
264
265
266
267
268
    }

    /**
     * sigma ratio between the block area in the x, y plane and D2
     */
    private get sigma(): number {
269
        if (this._cache.sigma === undefined) {
270
            if (this.prms.Cd0.v >= 2) {
271
272
273
274
                this._cache.sigma = 1;
            } else {
                this._cache.sigma = Math.PI / 4;
            }
Dorchies David's avatar
Dorchies David committed
275
        }
276
        return this._cache.sigma;
Dorchies David's avatar
Dorchies David committed
277
278
279
    }

    private get R(): number {
280
281
        if (this._cache.R === undefined) {
            this._cache.R = (1 - this.sigma * this.prms.C.v);
Dorchies David's avatar
Dorchies David committed
282
        }
283
        return this._cache.R;
Dorchies David's avatar
Dorchies David committed
284
285
286
287
    }

    /**
     * Bed friction coefficient Equation (3) (Cassan et al., 2016)
288
     * @param Y Water depth (m)
Dorchies David's avatar
Dorchies David committed
289
     */
290
    private calcCf(Y: number): number {
Dorchies David's avatar
Dorchies David committed
291

292
        if (this.prms.Ks.v < 1E-9) {
293
            // Between Eq (8) and (9) (Cassan et al., 2016)
294
            const reynolds = this.U0 * this.prms.Y.v / MacroRugo.nu;
295
            return 0.3164 / 4. * Math.pow(reynolds, -0.25);
Dorchies David's avatar
Dorchies David committed
296
297
        } else {
            // Equation (3) (Cassan et al., 2016)
298
            return 2 / Math.pow(5.1 * Math.log10(Y / this.prms.Ks.v) + 6, 2);
Dorchies David's avatar
Dorchies David committed
299
300
301
        }
    }

302
303
304
305
306
307
308
309
310
311
312
    /**
     * Interpolation of Cd0 for Cd from calibrated Cd0 (See #291)
     * Cx = 1.1 for Cd0 = 1 and Cx = 2.592 for Cd0 = 2
     */
    private get Cx(): number {
        if (this._cache.Cx === undefined) {
            this._cache.Cx = this.prms.Cd0.v * 1.4917 -0.3914;
        }
        return this._cache.Cx;
    }

Dorchies David's avatar
Dorchies David committed
313
314
315
    /**
     * Calculation of Cd : drag coefficient of a block under the actual flow conditions
     */
316
    private get Cd(): number {
Dorchies David's avatar
Dorchies David committed
317
        if (this._cache.Cd === undefined) {
318
            this._cache.Cd = this.Cx * Math.min(3, (1 + 1 / Math.pow(this.prms.Y.v / this.prms.PBD.v, 2)));
Dorchies David's avatar
Dorchies David committed
319
320
        }
        return this._cache.Cd;
Dorchies David's avatar
Dorchies David committed
321
322
323
324
325
326
327
    }

    /**
     * Calcul de Beta force ratio between drag and turbulent stress (Cassan et al. 2016 eq(8))
     * \Beta = (k / alpha_t) (C_d C k / D) / (1 - \sigma C)
     * @param alpha \alpha_t turbulent length scale (m) within the blocks layer
     */
328
    private calcBeta(alpha: number): number {
Dorchies David's avatar
Dorchies David committed
329
        return Math.min(100, Math.sqrt(this.prms.PBH.v * this.CdChD / alpha / this.R));
Dorchies David's avatar
Dorchies David committed
330
331
332
333
334
335
336
337
338
339
    }

    /**
     * Averaged velocity at a given vertical position (m.s-1)
     * @param alpha turbulent length scale (m) within the blocks layer
     * @param z dimensionless vertical position z / k
     */
    private calcUz(alpha: number, z: number = 1): number {
        const beta = this.calcBeta(alpha);
        // Equation (9) Cassan et al., 2016
340
341
342
        return this.u0 * Math.sqrt(
            beta * (this.prms.Y.v / this.prms.PBH.v - 1) * Math.sinh(beta * z) / Math.cosh(beta) + 1
        );
Dorchies David's avatar
Dorchies David committed
343
344
345
    }

    private get ustar(): number {
346
347
        if (this._cache.ustar === undefined) {
            this._cache.ustar = Math.sqrt(MacroRugo.g * this.prms.If.v * (this.prms.Y.v - this.prms.PBH.v));
Dorchies David's avatar
Dorchies David committed
348
        }
349
        return this._cache.ustar;
Dorchies David's avatar
Dorchies David committed
350
351
    }

352
    private resolveAlpha_t(alpha: number): number {
Dorchies David's avatar
Dorchies David committed
353
        /** s: minimum distance between blocks */
354
        const s = this.prms.PBD.v * (1 / Math.sqrt(this.prms.C.v) - 1);
Dorchies David's avatar
Dorchies David committed
355
356
357
358
359
360
        /** Equation(11) Cassan et al., 2016 */
        const l0 = Math.min(s, 0.15 * this.prms.PBH.v);
        // Equation(10) Cassan et al., 2016
        return alpha * this.calcUz(alpha) - l0 * this.ustar;
    }

361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
    private resolveQSubmerged(): number {
        /** Tirant d'eau (m) */
        const h: number = this.prms.Y.v;

        /** Concentration de blocs (-) */
        const C: number = this.prms.C.v;
        /** Paramètre de bloc : Diamètre (m) */
        const D: number = this.prms.PBD.v;
        /** Paramètre de bloc : Hauteur (m) */
        const k: number = this.prms.PBH.v;
        /** Slope (m/m) */
        const S: number = this.prms.If.v;
        /** Accélération gravité (m/s²) */
        const g = MacroRugo.g;
        /** Constante von Karman */
        const kappa = 0.41;
        /** Velocity at the bed §2.3.2 Cassan et al., 2016 */
        this.u0 = Math.sqrt(k * 2 * g * S * this.R
            / (this.Cd * C * k / D + this.calcCf(k) * this.R));
        /** turbulent length scale (m) within the blocks layer (alpha_t) */
        const alpha = uniroot(this.resolveAlpha_t, this, 1E-3, 10);
        /** averaged velocity at the top of blocks (m.s-1) */
        const uk = this.calcUz(alpha);
        /** Equation (13) Cassan et al., 2016 */
        const d = k - 1 / kappa * alpha * uk / this.ustar;
        /** Equation (14) Cassan et al., 2016 */
        const z0 = (k - d) * Math.exp(- kappa * uk / this.ustar);
        /** Integral of Equation (12) Cassan et al., 2016 */
        // tslint:disable-next-line:variable-name
        let Qsup: number;
        if (z0 > 0) {
            Qsup = this.ustar / kappa * (
                (h - d) * (Math.log((h - d) / z0) - 1)
                - ((k - d) * (Math.log((k - d) / z0) - 1))
            );
        } else {
            Qsup = 0;
        }

        // calcul intégrale dans la canopée----
        // tslint:disable-next-line:variable-name
        let Qinf: number = this.u0;
        let u = 0;
        let uOld: number;
        const step = 0.01;
        const zMax = 1 + step / 2;
        for (let z = step; z < zMax; z += step) {
            uOld = u;
            u = this.calcUz(alpha, z);
            Qinf += (uOld + u);
        }
        Qinf = Qinf / 2 * step * k;

        // Calcul de u moyen
        return this.U0 - (Qinf + Qsup) / h;
    }

Dorchies David's avatar
Dorchies David committed
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
    private resolveQEmergent(): number {
        // tslint:disable-next-line: variable-name
        const Cd = this.Cd * this.fFr;
        /** N from Cassan 2016 eq(2) et Cassan 2014 eq(12) */
        const N = (1 * this.calcCf(this.prms.Y.v)) / (this.prms.Y.v / this.prms.PBD.v * Cd * this.prms.C.v);

        // const U0i = Math.sqrt(
        //     2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v *
        //     (1 - this.sigma * this.prms.C.v) / (Cd * this.prms.C.v * (1 + N))
        // );

        return this.U0 - uniroot(this.resolveU0Complete, this, 1E-6, 100);
    }

    private resolveU0Complete(U0: number): number {

        const fFr = this.CalcfFr(U0);

436
        const alpha = 1 - Math.pow(1 * this.prms.C.v, 0.5) - 0.5 * this.sigma * this.prms.C.v;
Dorchies David's avatar
Dorchies David committed
437

Dorchies David's avatar
Dorchies David committed
438
439
        const N = (alpha * this.calcCf(this.prms.Y.v)) /
            (this.prms.Y.v / this.prms.PBD.v * this.Cd * fFr * this.prms.C.v);
Dorchies David's avatar
Dorchies David committed
440

Dorchies David's avatar
Dorchies David committed
441
        return U0 - Math.sqrt(
442
            2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v *
Dorchies David's avatar
Dorchies David committed
443
            (1 - this.sigma * this.prms.C.v) / (this.Cd * fFr * this.prms.C.v * (1 + N))
444
        );
Dorchies David's avatar
Dorchies David committed
445
446
    }

447
    /**
448
     * Calcul du ratio entre la vitesse moyenne à l'aval d'un block et la vitesse maximale
449
     * r = 1.1 pour un plot circulaire Cd0​=1.1 et r = 1.5 pour un plot à face plane Cd0​=2.6
450
     * Voir #283
451
452
453
     */
    private get r(): number {
        if (this._cache.r === undefined) {
454
            this._cache.r = 0.4 * this.prms.Cd0.v + 0.7;
455
456
457
458
459
460
461
        }
        return this._cache.r;
    }

    /**
     * Froude correction function (Cassan et al. 2014, Eq. 19)
     */
Dorchies David's avatar
Dorchies David committed
462
463
464
465
466
467
468
    private get fFr(): number {
        if (this._cache.fFr === undefined) {
            this._cache.fFr = this.CalcfFr(this.U0);
        }
        return this._cache.fFr;
    }

Dorchies David's avatar
Dorchies David committed
469
    /**
470
     * Calculation of Froude correction function (Cassan et al. 2014, Eq. 19)
Dorchies David's avatar
Dorchies David committed
471
     */
Dorchies David's avatar
Dorchies David committed
472
473
474
475
476
    private CalcfFr(U0: number): number {
        // tslint:disable-next-line:variable-name
        const Fr = U0 /
            (1 - Math.sqrt(MacroRugo.fracAxAy * this.prms.C.v)) /
            Math.sqrt(MacroRugo.g * this.prms.Y.v);
Dorchies David's avatar
Dorchies David committed
477

Dorchies David's avatar
Dorchies David committed
478
        if (Fr < 1.3) {
479
            return Math.pow(Math.min(1 / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2);
Dorchies David's avatar
Dorchies David committed
480
481
        } else {
            return Math.pow(Fr, -4 / 3);
Dorchies David's avatar
Dorchies David committed
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
        }
    }
}

/**
 * Searches the interval from <tt>lowerLimit</tt> to <tt>upperLimit</tt>
 * for a root (i.e., zero) of the function <tt>func</tt> with respect to
 * its first argument using Brent's method root-finding algorithm.
 *
 * Translated from zeroin.c in http://www.netlib.org/c/brent.shar.
 *
 * Copyright (c) 2012 Borgar Thorsteinsson <borgar@borgar.net>
 * MIT License, http://www.opensource.org/licenses/mit-license.php
 *
 * @param {function} func function for which the root is sought.
 * @param {number} lowerlimit the lower point of the interval to be searched.
 * @param {number} upperlimit the upper point of the interval to be searched.
 * @param {number} errorTol the desired accuracy (convergence tolerance).
 * @param {number} maxIter the maximum number of iterations.
 * @returns an estimate for the root within accuracy.
 *
 */
504
function uniroot<T>(func: (param: number) => number, thisArg: T, lowerLimit: number, upperLimit: number,
Mathias Chouet's avatar
Mathias Chouet committed
505
                    errorTol: number = 0, maxIter: number = 1000
Dorchies David's avatar
Dorchies David committed
506
) {
507
508
509
510
511
512
513
514
515
    let a = lowerLimit;
    let b = upperLimit;
    let c = a;
    let fa = func.call(thisArg, a);
    let fb = func.call(thisArg, b);
    let fc = fa;
    let tolAct;   // Actual tolerance
    let newStep;  // Step at this iteration
    let prevStep; // Distance from the last but one to the last approximation
516
    let p;        // Interpolation step is calculated in the form p/q; division is delayed until the last moment
517
    let q;
Dorchies David's avatar
Dorchies David committed
518
519
520

    while (maxIter-- > 0) {

521
        prevStep = b - a;
Dorchies David's avatar
Dorchies David committed
522
523
524
525
526
527
528

        if (Math.abs(fc) < Math.abs(fb)) {
            // Swap data for b to be the best approximation
            a = b, b = c, c = a;
            fa = fb, fb = fc, fc = fa;
        }

529
530
        tolAct = 1e-15 * Math.abs(b) + errorTol / 2;
        newStep = (c - b) / 2;
Dorchies David's avatar
Dorchies David committed
531

532
        if (Math.abs(newStep) <= tolAct || fb === 0) {
Dorchies David's avatar
Dorchies David committed
533
534
535
536
            return b; // Acceptable approx. is found
        }

        // Decide if the interpolation can be tried
537
        if (Math.abs(prevStep) >= tolAct && Math.abs(fa) > Math.abs(fb)) {
Dorchies David's avatar
Dorchies David committed
538
            // If prev_step was large enough and was in true direction, Interpolatiom may be tried
539
540
541
            let t1;
            let cb;
            let t2;
Dorchies David's avatar
Dorchies David committed
542
543
544
545
546
            cb = c - b;
            if (a === c) { // If we have only two distinct points linear interpolation can only be applied
                t1 = fb / fa;
                p = cb * t1;
                q = 1.0 - t1;
547
            } else { // Quadric inverse interpolation
Dorchies David's avatar
Dorchies David committed
548
549
550
551
552
553
554
                q = fa / fc, t1 = fb / fc, t2 = fb / fa;
                p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1));
                q = (q - 1) * (t1 - 1) * (t2 - 1);
            }

            if (p > 0) {
                q = -q;  // p was calculated with the opposite sign; make p positive
555
            } else {
Dorchies David's avatar
Dorchies David committed
556
557
558
                p = -p;  // and assign possible minus to q
            }

559
560
            if (p < (0.75 * cb * q - Math.abs(tolAct * q) / 2) &&
                p < Math.abs(prevStep * q / 2)) {
Dorchies David's avatar
Dorchies David committed
561
                // If (b + p / q) falls in [b,c] and isn't too large it is accepted
562
                newStep = p / q;
Dorchies David's avatar
Dorchies David committed
563
564
565
566
567
            }

            // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent
        }

568
569
        if (Math.abs(newStep) < tolAct) { // Adjust the step to be not less than tolerance
            newStep = (newStep > 0) ? tolAct : -tolAct;
Dorchies David's avatar
Dorchies David committed
570
571
572
        }

        a = b, fa = fb;     // Save the previous approx.
573
        b += newStep, fb = func.call(thisArg, b);  // Do step to a new approxim.
Dorchies David's avatar
Dorchies David committed
574
575
576
577
578
579
580

        if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) {
            c = a, fc = fa; // Adjust c for it to have a sign opposite to that of b
        }
    }
    return undefined;

581
}