p.[k] <- p.[k] / p.[0]
p.[0] <- 1.0
let e = 0.25 * p.[1]
- let b = ref (2.0 * e)
- let c = ref (!b ** 2.0)
- let mutable d = 0.75 * !c
- b := p.[3] + !b *(!c - p.[2])
+ let mutable b = 2.0 * e
+ let mutable c = b ** 2.0
+ let mutable d = 0.75 * c
+ b <- p.[3] + b *(c - p.[2])
let mutable a = p.[2] - d
- c := p.[4] + e * (e * a - p.[3])
+ c <- p.[4] + e * (e * a - p.[3])
a <- a - d
- let quadExecuted = ref false
+ let mutable quadExecuted = false
let quad () =
- if not !quadExecuted
+ if not quadExecuted
then
- p.[2] <- !c / !b
+ p.[2] <- c / b
quadroots p r
for k in 1..2 do
for j in 1..2 do
r.[j, k+2] <- r.[j, k]
p.[1] <- -p.[1]
- p.[2] <- !b
+ p.[2] <- b
quadroots p r
for k in 1..4 do
r.[1,k] <- r.[1,k] - e
- quadExecuted := true
+ quadExecuted <- true
p.[1] <- 0.5 * a
- p.[2] <- (p.[1] * p.[1] - !c) * 0.25
- p.[3] <- !b * !b / -64.0
+ p.[2] <- (p.[1] * p.[1] - c) * 0.25
+ p.[3] <- b * b / -64.0
if p.[3] < 0.0
then
cubicroots p r
then
d <- r.[1, k] * 4.0
a <- a + d
- if a >= 0.0 && !b >= 0.0
+ if a >= 0.0 && b >= 0.0
then
p.[1] <- sqrt d
- elif a <= 0.0 && !b <= 0.0
+ elif a <= 0.0 && b <= 0.0
then
p.[1] <- sqrt d
else
p.[1] <- -(sqrt d)
- b := 0.5 * (a + !b / p.[1])
+ b <- 0.5 * (a + b / p.[1])
quad ()
k <- 4
k <- k + 1
- if not !quadExecuted && p.[2] < 0.0
+ if not quadExecuted && p.[2] < 0.0
then
- b := sqrt !c
- d <- !b + !b - a
+ b <- sqrt c
+ d <- b + b - a
p.[1] <- 0.0
if d > 0.0
then
p.[1] <- sqrt d
- elif not !quadExecuted
+ elif not quadExecuted
then
if p.[1] > 0.0
then
- b := (sqrt p.[2]) * 2.0 + p.[1]
+ b <- (sqrt p.[2]) * 2.0 + p.[1]
else
- b := -(sqrt p.[2]) * 2.0 + p.[1]
+ b <- -(sqrt p.[2]) * 2.0 + p.[1]
- if !b <> 0.0
+ if b <> 0.0
then
p.[1] <- 0.0
else
for k in 1..4 do
r.[1, k] <- -e
r.[2, k] <- 0.0
- quadExecuted := true
+ quadExecuted <- true
quad ()
// Return a tuple (area, x intersections, y intersections)
-let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * float[]) option =
- let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha
- let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha
+let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option =
+ let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha
+ let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha
if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
then
printf "nroots = %d\n" nroots
#endif
- let ychk = [|
- for i in 1 .. nroots do
- if abs r.[2, i] < EPS
- then
- yield r.[1, i] * b1
+ let ychk = Array.init nroots (fun _ -> Double.MaxValue)
+ let mutable nychk = 0
+ for i in 1 .. nroots do
+ if abs r.[2, i] < EPS
+ then
+ ychk.[nychk] <- r.[1, i] * b1
+ nychk <- nychk + 1
#if DEBUG_LOG
- printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
+ printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
#endif
- |]
Array.sortInPlace ychk
#if DEBUG_LOG
printf "\t j=%d, ychk=%f\n" j ychk.[j]
#endif
- let nychk = Array.length ychk
let mutable nintpts = 0
let xint = Array.zeroCreate 4
| 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| _ -> -1.0
if nintpts = 0
- then Some (area, [||], [||])
+ then Some (float32 area, [||], [||])
else
- let xTransform = Array.zeroCreate nintpts
- let yTransform = Array.zeroCreate nintpts
+ let xTransform : float32[] = Array.zeroCreate nintpts
+ let yTransform : float32[] = Array.zeroCreate nintpts
for i in 0 .. (nintpts - 1) do
- xTransform.[i] <- cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
- yTransform.[i] <- sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
- Some (area, xTransform, yTransform)
\ No newline at end of file
+ xTransform.[i] <- float32 <| cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
+ yTransform.[i] <- float32 <| sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
+ Some (float32 area, xTransform, yTransform)
\ No newline at end of file