5 let private EPS = 1.0e-5
7 let inline private ellipse2tr (x
: float) (y
: float) (aa
: float) (bb
: float) (cc
: float) (dd
: float) (ee
: float) (ff
: float) : float =
8 aa
* x
* x
+ bb
* x
* y
+ cc
* y
* y
+ dd
* x
+ ee
* y
+ ff
10 let private nointpts
(a1
: float) (b1
: float) (a2
: float) (b2
: float) (h1
: float) (k1
: float) (h2
: float) (k2
: float) (phi_1
: float) (phi_2
: float) (h2_tr
: float) (k2_tr
: float) (aa
: float) (bb
: float) (cc
: float) (dd
: float) (ee
: float) (ff
: float) =
13 let area_1 = Math.PI * a1b1
14 let area_2 = Math.PI * a2b2
15 let relsize = a1b1 - a2b2
19 if (h2_tr
* h2_tr
) / (a1
* a1
) + (k2_tr
* k2_tr
) / (b1
* b1
) < 1.0
30 if abs
(h1
- h2
) < EPS && abs
(k1
- k2
) < EPS && abs
(area_1 - area_2) < EPS
34 type private PointType = TANGENT_POINT | INTERSECTION_POINT
36 let private istanpt
(x
: float) (y
: float) (a1
: float) (b1
: float) (aa
: float) (bb
: float) (cc
: float) (dd
: float) (ee
: float) (ff
: float) : PointType =
40 if x < 0.0 then -a1
else a1
45 then 2.0 * Math.PI - acos
(x / a1)
50 let x1 = a1 * cos
(theta + eps_radian)
51 let y1 = b1
* sin
(theta + eps_radian)
52 let x2 = a1 * cos
(theta - eps_radian)
53 let y2 = b1
* sin
(theta - eps_radian)
55 let test1 = ellipse2tr x1 y1 aa bb cc dd ee
ff
56 let test2 = ellipse2tr x2 y2 aa bb cc dd ee
ff
58 if test1 * test2 > 0.0
60 else INTERSECTION_POINT
63 let private twointpts
(x: float[]) (y
: float[]) (a1: float) (b1
: float) (phi_1
: float) (a2
: float) (b2
: float) (h2_tr
: float) (k2_tr
: float) (phi_2
: float) (aa
: float) (bb
: float) (cc
: float) (dd
: float) (ee
: float) (ff: float) =
65 then x.[0] <- if x.[0] < 0.0 then -a1 else a1
69 then 2.0 * Math.PI - acos (x.[0] / a1)
70 else acos (x.[0] / a1)
73 then x.[1] <- if x.[1] < 0.0 then -a1 else a1
77 then 2.0 * Math.PI - acos (x.[1] / a1)
78 else acos (x.[1] / a1)
86 let xmid = a1 * cos
((theta1 + theta2) / 2.0)
87 let ymid = b1
* sin
((theta1 + theta2) / 2.0)
89 if ellipse2tr xmid ymid aa bb cc dd ee
ff > 0.0
97 theta1 <- theta1 - 2.0 * Math.PI
99 let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
101 let mutable area1 = 0.5 * (a1 * b1
* (theta2 - theta1) + trsign * abs
(x.[0] * y
.[1] - x.[1] * y
.[0]))
105 area1 <- area1 + a1 * b1
107 let cosphi = cos
(phi_1
- phi_2
)
108 let sinphi = sin
(phi_1
- phi_2
)
110 let mutable x1_tr = (x.[0] - h2_tr
) * cosphi + (y
.[0] - k2_tr
) * -sinphi
111 let mutable y1_tr = (x.[0] - h2_tr
) * sinphi + (y
.[0] - k2_tr
) * cosphi
112 let mutable x2_tr = (x.[1] - h2_tr
) * cosphi + (y
.[1] - k2_tr
) * -sinphi
113 let mutable y2_tr = (x.[1] - h2_tr
) * sinphi + (y
.[1] - k2_tr
) * cosphi
117 x1_tr <- if x1_tr < 0.0 then -a2
else a2
121 theta1 <- 2.0 * Math.PI - acos (x1_tr / a2)
123 theta1 <- acos (x1_tr / a2)
127 x2_tr <- if x2_tr < 0.0 then -a2 else a2
131 theta2 <- 2.0 * Math.PI - acos (x2_tr / a2)
133 theta2 <- acos (x2_tr / a2)
141 let xmid = a2 * cos
((theta1 + theta2) / 2.0)
142 let ymid = b2
* sin
((theta1 + theta2) / 2.0)
144 let cosphi = cos
(phi_2
- phi_1
)
145 let sinphi = sin
(phi_2
- phi_1
)
146 let xmid_rt = xmid * cosphi + ymid * -sinphi + h2_tr
147 let ymid_rt = xmid * sinphi + ymid * cosphi + k2_tr
149 if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1
* b1
) > 1.0
157 theta1 <- theta1 - 2.0 * Math.PI
159 let trsign = if theta2 - theta1 > Math.PI then 1.0 else -1.0
161 let mutable area2 = 0.5 * (a2 * b2
* (theta2 - theta1) + trsign * abs
(x1_tr * y2_tr - x2_tr * y1_tr))
164 area2 <- area2 + a2 * b2
169 let private threeintpts
(xint
: float[]) (yint
: float[]) (a1: float) (b1
: float) (phi_1
: float) (a2: float) (b2
: float) (h2_tr
: float) (k2_tr
: float) (phi_2
: float) (aa
: float) (bb
: float) (cc
: float) (dd
: float) (ee
: float) (ff: float) : float =
170 let mutable tanpts = 0
171 let mutable tanindex = 0
173 if istanpt
xint.[i
] yint
.[i
] a1 b2 aa bb cc dd ee
ff = TANGENT_POINT
175 tanpts <- tanpts + 1;
191 twointpts
xint yint
a1 b1 phi_1
a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee
ff
194 let private fourintpts
(xint: float[]) (yint
: float[]) (a1: float) (b1
: float) (phi_1
: float) (a2: float) (b2
: float) (h2_tr
: float) (k2_tr
: float) (phi_2
: float) (aa
: float) (bb
: float) (cc
: float) (dd
: float) (ee
: float) (ff: float) : float =
197 let area_1 = Math.PI * a1b1
198 let area_2 = Math.PI * a2b2
200 let theta = Array.zeroCreate
4
205 xint.[i
] <- if xint.[i
] < 0.0 then -a1 else a1
206 theta.[i
] <- if yint
.[i
] < 0.0 then 2.0 * Math.PI - acos (xint.[i
] / a1) else acos (xint.[i
] / a1)
213 let mutable k = j
- 1
219 theta.[k+1] <- theta.[k]
220 xint.[k+1] <- xint.[k]
221 yint
.[k+1] <- yint
.[k]
228 let area1 = 0.5 * abs
((xint.[2] - xint.[0]) * (yint
.[3] - yint
.[1]) - (xint.[3] - xint.[1]) * (yint
.[2] - yint
.[0]))
230 let cosphi = cos
(phi_1
- phi_2
)
231 let sinphi = sin
(phi_1
- phi_2
)
233 let theta_tr = Array.zeroCreate
4
234 let xint_tr = Array.zeroCreate
4
235 let yint_tr = Array.zeroCreate
4
238 xint_tr.[i
] <- (xint.[i
] - h2_tr
) * cosphi + (yint
.[i
] - k2_tr
) * -sinphi
239 yint_tr.[i
] <- (xint.[i
] - h2_tr
) * sinphi + (yint
.[i
] - k2_tr
) * cosphi
241 if abs
xint_tr.[i
] > a2
243 xint_tr.[i
] <- if xint_tr.[i
] < 0.0 then -a2 else a2
245 theta_tr.[i
] <- if yint_tr.[i
] < 0.0 then 2.0 * Math.PI - acos (xint_tr.[i
] / a2) else acos (xint_tr.[i
] / a2)
247 let xmid = a1 * cos
((theta.[0] + theta.[1]) / 2.0)
248 let ymid = b1
* sin
((theta.[0] + theta.[1]) / 2.0)
250 let mutable area2, area3
, area4
, area5
= 0.0, 0.0, 0.0, 0.0
252 if ellipse2tr xmid ymid aa bb cc dd ee
ff < 0.0
254 area2 <- 0.5 * (a1b1 * (theta.[1] - theta.[0]) - abs
(xint.[0] * yint
.[1] - xint.[1] * yint
.[0]))
255 area4
<- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs
(xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1]))
257 if theta_tr.[3] > theta_tr.[0]
259 area5
<- 0.5 * (a2b2 * (theta_tr.[0] - (theta_tr.[3] - 2.0 * Math.PI)) - abs
(xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
261 area5
<- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs
(xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
263 area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs
(xint.[1] * yint
.[2] - xint.[2] * yint
.[1]))
264 area3
<- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs
(xint.[3] * yint
.[0] - xint.[0] * yint
.[3]))
265 area4
<- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs
(xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
266 area5
<- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs
(xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))
270 area5
<- area5
+ area_2
274 area4
<- area4
+ area_2
278 area3
<- area3
+ area_1
282 area2 <- area2 + area_1
284 area1 + area2 + area3
+ area4
+ area5
287 let private quadroots
(p
: float[]) (r
: float[,]) =
288 let mutable b = -p
.[1] / (2.0 * p
.[0])
289 let c = p
.[2] / p
.[0]
290 let mutable d = b * b - c
311 let private cubicroots
(p
: float[]) (r
: float[,]) =
314 p
.[k] <- p
.[k] / p
.[0]
317 let mutable t = s * p
.[1]
318 let mutable b = 0.5 * (s * (t / 1.5 - p
.[2]) + p
.[3])
319 t <- (t - p
.[2]) / 3.0
320 let mutable c = t * t * t
321 let mutable d = b * b - c
325 d <- ((sqrt
d) + (abs
b)) ** (1.0 / 3.0)
332 d <- sqrt
(0.75) * (b - c)
337 if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0
350 then d <- (atan
1.0) / 1.5
351 else d <- atan
((sqrt
-d) / (abs
b)) / 3.0
354 then b <- 2.0 * (sqrt
t)
355 else b <- -2.0 * (sqrt
t)
358 t <- -(sqrt
0.75
) * (sin
d) * b - 0.5 * c
382 let private biquadroots
(p
: float[]) (r
: float[,]) =
386 p
.[k] <- p
.[k] / p
.[0]
389 let b = ref (2.0 * e)
390 let c = ref (!b ** 2.0)
391 let mutable d = 0.75 * !c
392 b := p
.[3] + !b *(!c - p
.[2])
393 let mutable a = p
.[2] - d
394 c := p
.[4] + e * (e * a - p
.[3])
397 let quadExecuted = ref false
405 r
.[j
, k+2] <- r
.[j
, k]
410 r
.[1,k] <- r
.[1,k] - e
414 p.[2] <- (p.[1] * p.[1] - !c) * 0.25
415 p.[3] <- !b * !b / -64.0
421 if r
.[2, k] = 0.0 && r
.[1, k] > 0.0
425 if a >= 0.0 && !b >= 0.0
428 elif
a <= 0.0 && !b <= 0.0
433 b := 0.5 * (a + !b / p.[1])
438 if not
!quadExecuted && p.[2] < 0.0
446 elif not
!quadExecuted
450 b := (sqrt
p.[2]) * 2.0 + p.[1]
452 b := -(sqrt
p.[2]) * 2.0 + p.[1]
468 let EEOverlapArea (e1
: Ellipse) (e2
: Ellipse) : float =
469 let { cx
= h1
; cy
= k1
; a = a1; b = b1
; alpha
= phi_1
} = e1
470 let { cx
= h2
; cy
= k2
; a = a2; b = b2
; alpha
= phi_2
} = e2
472 if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
476 let phi_1 = phi_1 % Math.PI
477 let phi_2 = phi_2 % Math.PI
478 let h2_tr, k2_tr
, phi_2r
=
479 let cosphi = cos
phi_1
480 let sinphi = sin
phi_1
481 (h2
- h1
) * cosphi + (k2
- k1
) * sinphi, (h1
- h2
) * sinphi + (k2
- k1
) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
483 let cosphi = cos
phi_2r
484 let cosphi2 = cosphi ** 2.0
485 let sinphi = sin
phi_2r
486 let sinphi2 = sinphi ** 2.0
487 let cosphisinphi = 2.0 * cosphi * sinphi
490 let tmp0 = (cosphi * h2_tr + sinphi * k2_tr
) / a22
491 let tmp1 = (sinphi * h2_tr - cosphi * k2_tr
) / b22
492 let tmp2 = cosphi * h2_tr + sinphi * k2_tr
493 let tmp3 = sinphi * h2_tr - cosphi * k2_tr
495 let aa = cosphi2 / a22 + sinphi2 / b22
496 let bb = cosphisinphi / a22 - cosphisinphi / b22
497 let cc = sinphi2 / a22 + cosphi2 / b22
498 let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1
499 let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1
500 let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0
503 (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff)
504 2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff)
505 a1 * a1 * ((b1 * b1 * (2.0 * aa * cc - bb * bb) + dd * dd - 2.0 * aa * ff) - 2.0 * a1 * a1 * aa * aa) + b1 * b1 * (2.0 * cc * ff + ee * ee)
506 2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee))
507 a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
510 let py = Array.zeroCreate
<float> 5
511 let r = Array2D.zeroCreate
<float> 3 5
517 py.[4-i
] <- cy.[i
] / cy.[4]
522 elif abs
cy.[3] > EPS
525 py.[3 - i
] <- cy.[i
] / cy.[3]
530 elif abs
cy.[2] > EPS
533 py.[2-i
] <- cy.[i
] / cy.[2]
538 elif abs
cy.[1] > EPS
540 r.[1, 1] <- -cy.[0] / cy.[1]
548 for i
in 1 .. nroots do
549 if abs
r.[2, i
] < EPS
553 Array.sortInPlace
ychk
555 let nychk = Array.length
ychk
556 let mutable nintpts = 0
558 let xint = Array.zeroCreate
4
559 let yint = Array.zeroCreate
4
561 let mutable returnValue = 0.0
564 while returnValue = 0.0 && i < nychk do
565 if not
(i < nychk - 1 && abs
(ychk.[i] - ychk.[i+1]) < EPS / 2.0)
567 let x1 = if abs
ychk.[i] > b1 then 0.0 else a1 * sqrt
(1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
570 if abs
(ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
572 nintpts <- nintpts + 1
577 xint.[nintpts-1] <- x1
578 yint.[nintpts-1] <- ychk.[i]
580 if returnValue <> -1.0 && abs
(ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs
(x2 - x1) > EPS
582 nintpts <- nintpts + 1
587 xint.[nintpts-1] <- x2
588 yint.[nintpts-1] <- ychk.[i]
592 if returnValue = -1.0
597 | 0 | 1 -> nointpts
a1 b1 a2 b2 h1 k1 h2 k2
phi_1 phi_2 h2_tr k2_tr
aa bb cc dd ee ff
598 | 2 -> match istanpt
xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
599 | TANGENT_POINT -> nointpts
a1 b1 a2 b2 h1 k1 h2 k2
phi_1 phi_2 h2_tr k2_tr
aa bb cc dd ee ff
600 | INTERSECTION_POINT -> twointpts
xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr
phi_2 aa bb cc dd ee ff
601 | 3 -> threeintpts
xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr
phi_2 aa bb cc dd ee ff
602 | 4 -> fourintpts
xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr
phi_2 aa bb cc dd ee ff