Commit 6e4eb30a authored by Dorchies David's avatar Dorchies David
Browse files

Merge branch 'devel'

parents 2814d644 f0a8431a
function [f] = calcfFr(Frg)
if abs(Frg)>1.3
f=(Frg.^(-2/3)).^2;
function [f] = calcfFr(Frg, rQ)
if ~exists("rQ","local") then
rQ = 0.4 * Cd0 + 0.7
end
if abs(Frg) < 1
f=(min(rQ/(1-(Frg.^2)/4),abs(Frg).^(-2/3))).^2;
else
f=(min(1./(1-(Frg.^2)/4),abs(Frg).^(-2/3))).^2;
f = 1
end
endfunction
......@@ -19,7 +19,7 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
fFr = calcfFr(Frg)
//Cd=Cd0.*(1+0.4./(pf./D).^2);
Cd=Cd0.*min((1+1./(pf./D).^2), 3);
Cd=Cd0.*min(1E9,(1+1./(pf./D).^2));
//Cd=Cd0.*(0.8-2*C).*(1+0.4./(pf./D).^2).*fFr;
R=(1-sigma*C);//%.*(1-C.^0.5).^2;
if(bDbg) then
......@@ -88,8 +88,6 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
else
alpha=1-(1.*C).^0.5-1/2*sigma.*C;
if bDbg then print_r('alpha'); end;
Cd2=Cd.*fFr;
if bDbg then print_r('Cd2'); end;
......@@ -104,15 +102,17 @@ function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg)
end
if bDbg then print_r('Cf'); end;
N= (1 .*Cf)./(pf./D.*Cd2.*C);
alpha=1-(1.*C).^0.5-1/2*sigma.*C;
if bDbg then print_r('alpha'); end;
N= (alpha .*Cf)./(pf./D.*Cd2.*C);
if bDbg then print_r('N'); end;
U0a=(2*g.*S.*D.*(R)./(Cd2.*C.*(1+N))).^0.5;
if bDbg then print_r("U0a");end;
U0complet = find_U0_complet(U0a, pf,C,D,sigma,Cd0,Cf,S, bDbg)
U0complet = find_U0_complet(U0a, pf,C,D,sigma,Cd,Cf,S, bDbg)
if bDbg then print_r("U0complet");end;
[u fval] = fminsearch(list(find_U0_complet, pf,C,D,sigma,Cd0,Cf,S), U0a, opt)
[u fval] = fminsearch(list(find_U0_complet, pf,C,D,sigma,Cd,Cf,S), U0a, opt)
res=abs(U0-u);
......
function [res]=find_U0_complet(U01,h,C,D,sigma,Cd0,cfmean,S, bDbgU0)
function [res]=find_U0_complet(U01,h,C,D,sigma,Cd,cfmean,S, bDbgU0)
if ~exists("bDbgU0","local"); bDbgU0 = %f; end;
g=9.81;
Frexp=U01./(1-(1.*C).^0.5)./(9.81*h)^0.5;
//Frexp=U01./(9.81*h)^0.5;
fFr=(min(1./(1-(Frexp.^2)/4),Frexp.^(-2/3))).^2;
if Frexp>=1.5
fFr=(Frexp.^(-2/3)).^2;
end
fFr = calcfFr(Frexp)
alpha=1-(1.*C).^0.5-1/2*sigma.*C;
//Cd=Cd0.*(0.8-2*C).*(1+0.4./(h./D).^2).*fFr;
Cd=Cd0.*(1+1./(h./D).^2).*fFr;
Cd2=Cd.*fFr;
N= (alpha.*cfmean)./(h./D.*Cd.*C);
N= (alpha.*cfmean)./(h./D.*Cd2.*C);
if bDbgU0 then
print_r("Frexp")
......@@ -25,7 +19,7 @@ if bDbgU0 then
print_r("N")
end
res=abs(U01-(2*g.*S.*D.*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5);
res=abs(U01-(2*g.*S.*D.*(1-(sigma*C))./(Cd2.*C.*(1+N))).^0.5);
if bDbgU0 then
print_r("res")
end
......
......@@ -35,12 +35,12 @@ function macrorugo_resultComp(z_amont, S, long, Q, L, pf, C, Cd0, h, D)
q_technique=0.815*(pf/D)^1.45*S^0.557*C^(-0.456)*(9.81*D)^0.5.*D*L;
V_technique=4.54*(pf/D)^0.32*S^0.56*(9.81*D)^0.5;
end
if Cd0 > 2
coeff_contraction=1.5;
if Cd0 > 1.99
rv=1.5;
else
coeff_contraction=1.1;
rv=1.1;
end
Vmax=Vg.*coeff_contraction*calcfFr(Fr);
Vmax=Vg.*sqrt(calcfFr(Fr, rv));
print_r("Vmax");
print_r("V_technique");
end
......
......@@ -16,34 +16,33 @@ z_amont = 12.5; // Cote amont (m)
long = 6; // Longueur rampe (m)
// *****************************************************************************
printf("\n*** Emergent conditions Cd=1.1003***\n")
printf("\n*** Emergent conditions Cd=1.***\n")
// *****************************************************************************
h = 0.6
k = 0.7
Cd0 = 1.1003 // Forme ronde
h = 0.4
k = 0.6
Cd0 = 1. // Forme ronde
macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg)
// *****************************************************************************
printf("\n*** Emergent conditions Cd=2.592***\n")
printf("\n*** Emergent conditions Cd=2***\n")
// *****************************************************************************
h = 0.6
k = 0.7
Cd0 = 2.592 // Forme plane
h = 0.4
k = 0.6
Cd0 = 2. // Forme plane
macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg)
// *****************************************************************************
printf("\n*** Submerged conditions Cd=1.1003***\n")
printf("\n*** Submerged conditions Cd=1***\n")
// *****************************************************************************
k = 0.7
k = 0.6
h = 0.8
Cd0 = 1.1003 // Forme ronde
C = 0.13
Cd0 = 1. // Forme ronde
macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg)
// *****************************************************************************
printf("\n*** JalHyd #85 ***\n")
// *****************************************************************************
Cd0 = 1.1003
Cd0 = 1.
k = 0.8
C = 0.2; // Concentration
Q = []
......
......@@ -4,29 +4,29 @@ import { MacrorugoParams } from "../../src/macrorugo/macrorugo_params";
import { checkResult } from "../test_func";
/*
** Emergent conditions Cd=1.1003***
*** Emergent conditions Cd=1.***
*** INPUT ***
ks: 0.010000,
D: 0.500000,
k = 0.700000
Cd0: 1.100300,
k = 0.600000
Cd0: 1.000000,
S: 0.050000,
B: 1.000000,
h: 0.600000,
h: 0.400000,
C: 0.130000,
RESULTS:
find_Q_nat(1.052827)=0.600769
Q=0.692366 fVal=0.000000
find_Q_nat(0.736243)=0.913414
Q=0.370878 fVal=0.000000
ZF2: 12.200000,
Vdeb: 1.153943,
Vg: 1.804601,
Fr: 0.743826,
PV: 566.009124,
Vdeb: 0.927194,
Vg: 1.449999,
Fr: 0.731987,
PV: 454.788830,
flowcond: emergent,
Vmax: 2.673506,
V_technique: 1.991299,
q_technique: 0.561860,
Strickler: 7.254351,
Vmax: 1.785249,
V_technique: 1.748990,
q_technique: 0.312101,
Strickler: 7.637991,
*/
function macroRugoInstanceEmergentCd1(): MacroRugo {
......@@ -36,12 +36,12 @@ function macroRugoInstanceEmergentCd1(): MacroRugo {
6, // L
1, // B
0.05, // If
0.692366, // Q
0.6, // h
0.370878, // Q
0.4, // h
0.01, // Ks
0.13, // C
0.5, // D
0.7, // k
0.6, // k
1 // Cd0
)
);
......@@ -50,95 +50,95 @@ function macroRugoInstanceEmergentCd1(): MacroRugo {
const macroRugoExtraResultEmergentCd1: { [key: string]: number | MacroRugoFlowType } = {
ENUM_MacroRugoFlowType: MacroRugoFlowType.EMERGENT,
ZF2: 12.200000,
Vdeb: 1.153943,
Vg: 1.804601,
Fr: 0.743826,
PV: 566.009124,
Vmax: 2.673506,
Strickler: 7.254351
Vdeb: 0.927194,
Vg: 1.449999,
Fr: 0.731987,
PV: 454.788830,
Vmax: 1.785249,
Strickler: 7.637991
};
/*
** Emergent conditions Cd=2.592***
*** Emergent conditions Cd=2***
*** INPUT ***
ks: 0.010000,
D: 0.500000,
k = 0.700000
Cd0: 2.592000,
k = 0.600000
Cd0: 2.000000,
S: 0.050000,
B: 1.000000,
h: 0.600000,
h: 0.400000,
C: 0.130000,
RESULTS:
find_Q_nat(0.675214)=0.321449
Q=0.482344 fVal=0.000000
find_Q_nat(0.512451)=0.770391
Q=0.204295 fVal=0.000000
ZF2: 12.200000,
Vdeb: 0.803907,
Vg: 1.257195,
Fr: 0.518194,
PV: 394.316317,
Vdeb: 0.510737,
Vg: 0.798720,
Fr: 0.403209,
PV: 250.516618,
flowcond: emergent,
Vmax: 2.166970,
V_technique: 1.991299,
q_technique: 0.561860,
Strickler: 5.053822
*/
Vmax: 1.248838,
V_technique: 1.427752,
q_technique: 0.266857,
Strickler: 4.207323,
*/
function macroRugoInstanceEmergentCd2(): MacroRugo {
const nub: MacroRugo = macroRugoInstanceEmergentCd1();
nub.prms.Cd0.singleValue = 2;
nub.prms.Q.singleValue = 0.481698;
nub.prms.Q.singleValue = 0.204295;
return nub;
}
const macroRugoExtraResultEmergentCd2: { [key: string]: number | MacroRugoFlowType } = {
ENUM_MacroRugoFlowType: MacroRugoFlowType.EMERGENT,
ZF2: 12.200000,
Vdeb: 0.803907,
Vg: 1.257195,
Fr: 0.518194,
PV: 394.316317,
Vmax: 2.166970,
Strickler: 5.053822
Vdeb: 0.510737,
Vg: 0.798720,
Fr: 0.403209,
PV: 250.516618,
Vmax: 1.248838,
Strickler: 4.207323
};
/*
** Submerged conditions Cd=1.1003***
*** Submerged conditions Cd=1***
*** INPUT ***
ks: 0.010000,
D: 0.500000,
k = 0.700000
Cd0: 1.100300,
k = 0.600000
Cd0: 1.000000,
S: 0.050000,
B: 1.000000,
h: 0.800000,
C: 0.130000,
RESULTS:
find_Q_nat(1.403770)=0.143414
Q=1.289038 fVal=0.000000
find_Q_nat(1.472487)=0.095732
Q=1.549073 fVal=0.000000
ZF2: 12.200000,
Vdeb: 1.611298,
Vg: 2.519839,
Fr: 0.899484,
PV: 790.341606,
Vdeb: 1.936341,
Vg: 3.028159,
Fr: 1.080934,
PV: 949.775168,
flowcond: immerge,
q_technique: 0.940450,
Strickler: 8.361756
q_technique: 1.060933,
Strickler: 10.048552,
*/
function macroRugoInstanceSubmerged(): MacroRugo {
const nub: MacroRugo = macroRugoInstanceEmergentCd1();
nub.prms.Y.singleValue = 0.8;
nub.prms.Q.singleValue = 1.289038;
nub.prms.Q.singleValue = 1.549073;
return nub;
}
const macroRugoExtraResultSubmerged: { [key: string]: number | MacroRugoFlowType } = {
ENUM_MacroRugoFlowType: MacroRugoFlowType.SUBMERGED,
ZF2: 12.200000,
Vdeb: 1.611298,
PV: 790.341606,
Strickler: 8.361756
Vdeb: 1.936341,
PV: 949.775168,
Strickler: 10.048552,
};
function MacroRugoFactory(sInstance: string): MacroRugo {
......@@ -186,7 +186,7 @@ function testMacroRugoConfig(sInstance: string, Q0: number, fVal0: number, mrExt
// tslint:disable-next-line:no-string-literal
nubit["setFlowType"]();
// tslint:disable-next-line:no-string-literal
expect(nubit["resolveQ"]()).toBeCloseTo(fVal0, 5);
expect(Math.abs(nubit["resolveQ"]())).toBeCloseTo(fVal0, 5);
});
const nub = MacroRugoFactory(sInstance);
for (const prm of nub.prms) {
......@@ -215,8 +215,8 @@ function macroRugoInstanceJalHyd85(): MacroRugo {
describe("Class MacroRugo: ", () => {
testMacroRugoConfig("EmergentCd1", 1.052827, 0.600769, macroRugoExtraResultEmergentCd1);
testMacroRugoConfig("EmergentCd2", 0.675214, 0.321449, macroRugoExtraResultEmergentCd2);
testMacroRugoConfig("EmergentCd1", 0.736243, 0.913414, macroRugoExtraResultEmergentCd1);
testMacroRugoConfig("EmergentCd2", 0.512451, 0.770391, macroRugoExtraResultEmergentCd2);
// Le test passe en debug mais pas sous Jasmine !??
// describe(`Condition Submerged` , () => {
// it(`resolveAlpha_t(0.07) should be 0.074752`, () => {
......@@ -224,7 +224,7 @@ describe("Class MacroRugo: ", () => {
// expect(macroRugoInstanceSubmerged()["resolveAlpha_t"](0.07)).toBeCloseTo(0.074752, 5);
// });
// });
testMacroRugoConfig("Submerged", 1.403770, 0.143414, macroRugoExtraResultSubmerged);
testMacroRugoConfig("Submerged", 1.472487, 0.095732, macroRugoExtraResultSubmerged);
describe("JalHyd #85", () => {
it("CalcSerie Q should return the good result :)", () => {
......@@ -232,7 +232,7 @@ describe("Class MacroRugo: ", () => {
const nubMR2 = macroRugoInstanceJalHyd85();
nubMR.prms.Y.setValues(0.7, 1.2, 0.1);
const aQ: number[] = [0.688476, 0.824387, 1.166239, 1.475853, 1.850393, 2.288123];
const aQ: number[] = [0.662919, 0.793504, 1.221194, 1.542033, 1.928510, 2.378835];
nubMR.CalcSerie();
for (let i = 0; i < aQ.length; i++) {
......
......@@ -874,7 +874,8 @@ describe("vérificateur de franchissement −", () => {
mr.prms.B.singleValue = 1.789; // 2 patterns
mr.prms.C.singleValue = 0.2; // limite la vitesse max
mr.prms.PBH.singleValue = 0.6; // limite la vitesse max
mr.prms.Cd0.singleValue = 2.6; // limite la vitesse max
mr.prms.Cd0.singleValue = 2; // limite la vitesse max
mr.prms.If.singleValue = 0.04;
// vérificateur
const v = new Verificateur();
v.nubToVerify = mr;
......@@ -900,7 +901,7 @@ describe("vérificateur de franchissement −", () => {
// vérificateur
const v = new Verificateur();
v.nubToVerify = mr;
v.speciesList.push(FishSpecies[FishSpecies.SPECIES_3a]); // 2.533 > 2.50
v.speciesList.push(FishSpecies[FishSpecies.SPECIES_8a]); // > 1.50
// résultat
const res = v.CalcSerie();
expect(res.ok).toBe(false);
......
......@@ -39,6 +39,21 @@ export class MacroRugo extends FishPass {
private _cache: { [key: string]: number };
/** Coefficients used in f_h*(h*) */
private paramFhStar: [number, number, number] = [1, 1, 2];
/** Coefficient used in rQ */
private paramRQ: [number, number] = [0.4, 0.7];
/** Coefficient used in rQ */
private paramRV: [number, number] = [0.4, 0.7];
/** Maximum value for Cd */
private paramMaxCd: number = Infinity;
/** true: Cd0 * min(3, fh), false : min(6, Cd0 * fh) */
private paramCdNewVersion: boolean = true;
constructor(prms: MacrorugoParams, dbg: boolean = false) {
super(prms, dbg);
this._cache = {};
......@@ -117,7 +132,7 @@ export class MacroRugo extends FishPass {
}
r.resultElement.values.Fr = resFr;
// Vitesse maximale
r.resultElement.values.Vmax = this.r * r.resultElement.values.Vg * Math.sqrt(this.CalcfFr(resVdeb));
r.resultElement.values.Vmax = r.resultElement.values.Vg * Math.sqrt(this.CalcfFr(resVdeb, this.rV));
}
// Puissance dissipée
r.resultElement.values.PV = 1000 * MacroRugo.g * r.resultElement.values.Vdeb * this.prms.If.v;
......@@ -129,23 +144,6 @@ export class MacroRugo extends FishPass {
} else {
r.resultElement.values.ENUM_MacroRugoFlowType = MacroRugoFlowType.SUBMERGED;
}
// Vitesse et débit du guide technique
/* let cQ: [number, number, number, number];
let cV: [number, number, number];
let hdk: number;
if (r.resultElement.values.ENUM_MacroRugoFlowType === MacroRugoFlowType.SUBMERGED) {
cQ = [0.955, 2.282, 0.466, -0.23];
hdk = this.prms.PBH.v;
} else {
hdk = this.prms.PBD.v;
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];
cV = [4.54, 0.32, 0.56];
}
} */
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));
......@@ -158,7 +156,7 @@ export class MacroRugo extends FishPass {
public Equation(sVarCalc: string): Result {
const Q = this.prms.Q.v;
const q0 = Math.sqrt(2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v * (1 - (this.sigma * this.prms.C.v)) /
(this.Cx * this.prms.C.v)) * this.prms.Y.v * this.prms.B.v;
(this.prms.Cd0.v * this.prms.C.v)) * this.prms.Y.v * this.prms.B.v;
let r: Result;
if (q0 > 0) {
this.setFlowType();
......@@ -231,16 +229,13 @@ export class MacroRugo extends FishPass {
// Reset cached variables depending on Q (or not...)
this._cache = {};
/** adimensional water depth */
const hstar: number = this.prms.Y.v / this.prms.PBH.v;
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);
const a = (this.prms.Y.v / this.prms.PBH.v - 1) / (MacroRugo.limitSubmerg - 1);
return (1 - a) * this.resolveQEmergent() + a * this.resolveQSubmerged();
}
}
......@@ -299,23 +294,16 @@ export class MacroRugo extends FishPass {
}
}
/**
* 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;
}
/**
* Calculation of Cd : drag coefficient of a block under the actual flow conditions
*/
private get Cd(): number {
if (this._cache.Cd === undefined) {
this._cache.Cd = this.Cx * Math.min(3, (1 + 1 / Math.pow(this.prms.Y.v / this.prms.PBD.v, 2)));
if(this.paramCdNewVersion) {
this._cache.Cd = this.prms.Cd0.v * Math.min(this.paramMaxCd, (this.paramFhStar[0] + this.paramFhStar[1] / Math.pow(this.prms.Y.v / this.prms.PBD.v, this.paramFhStar[2])));
} else {
this._cache.Cd = Math.min(this.paramMaxCd, this.prms.Cd0.v * (this.paramFhStar[0] + this.paramFhStar[1] / Math.pow(this.prms.Y.v / this.prms.PBD.v, this.paramFhStar[2])));
}
}
return this._cache.Cd;
}
......@@ -416,44 +404,45 @@ export class MacroRugo extends FishPass {
}
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);
const alpha = 1 - Math.pow(1 * this.prms.C.v, 0.5) - 0.5 * this.sigma * this.prms.C.v;
// tslint:disable-next-line: variable-name
const Cd = this.Cd * this.CalcfFr(U0, this.rQ);
const N = (alpha * this.calcCf(this.prms.Y.v)) /
(this.prms.Y.v / this.prms.PBD.v * this.Cd * fFr * this.prms.C.v);
(this.prms.Y.v / this.prms.PBD.v * Cd * this.prms.C.v);
return U0 - Math.sqrt(
2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v *
(1 - this.sigma * this.prms.C.v) / (this.Cd * fFr * this.prms.C.v * (1 + N))
(1 - this.sigma * this.prms.C.v) / (Cd * this.prms.C.v * (1 + N))
);
}
/**
* Calcul du ratio entre la vitesse moyenne à l'aval d'un block et la vitesse maximale
* r = 1.1 pour un plot circulaire Cd0​=1.1 et r = 1.5 pour un plot à face plane Cd0​=2.6
* Voir #283
* r = 1.1 pour un plot circulaire Cd0​=1 et r = 1.5 pour un plot à face plane Cd0​=2
*/
private get r(): number {
if (this._cache.r === undefined) {
this._cache.r = 0.4 * this.prms.Cd0.v + 0.7;
private get rV(): number {
if (this._cache.rV === undefined) {
this._cache.rV = this.paramRV[0] * this.prms.Cd0.v + this.paramRV[1];
}
return this._cache.r;
return this._cache.rV;
}
/**
* Perte de charge supplémentaire due à la forme (voir fFr)
* r = 1 pour un plot circulaire Cd0​=1 et r = 1.25 pour un plot à face plane Cd0​=2
*/
private get rQ(): number {
if (this._cache.rQ === undefined) {
this._cache.rQ = this.paramRQ[0] * this.prms.Cd0.v + this.paramRQ[1];
}
return this._cache.rQ;
}
/**
......@@ -461,7 +450,7 @@ export class MacroRugo extends FishPass {
*/
private get fFr(): number {
if (this._cache.fFr === undefined) {
this._cache.fFr = this.CalcfFr(this.U0);
this._cache.fFr = this.CalcfFr(this.U0, this.rQ);
}
return this._cache.fFr;
}
......@@ -469,17 +458,12 @@ export class MacroRugo extends FishPass {
/**
* Calculation of Froude correction function (Cassan et al. 2014, Eq. 19)
*/
private CalcfFr(U0: number): number {
private CalcfFr(U0: number, r: 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);
if (Fr < 1.3) {
return Math.pow(Math.min(1 / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2);
} else {
return Math.pow(Fr, -4 / 3);
}
return Math.max(1, Math.pow(Math.min(r / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2));
}
}