let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff
let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff
+#if DEBUG_LOG
+ printf "\t\t--- debug istanpt with (x,y)=(%f, %f), A1=%f, B1=%f\n" x y a1 b1
+ printf "theta=%f\n" theta
+ printf "eps_Radian=%f\n" eps_radian
+ printf "(x1, y1)=(%f, %f)\n" x1 y1
+ printf "(x2, y2)=(%f, %f)\n" x2 y2
+ printf "test1=%f\n" test1
+ printf "test2=%f\n" test2
+#endif
+
if test1 * test2 > 0.0
then TANGENT_POINT
else INTERSECTION_POINT
if area1 < 0.0
then
+#if DEBUG_LOG
+ printf "TWO area1=%f\n" area1
+#endif
area1 <- area1 + a1 * b1
let cosphi = cos (phi_1 - phi_2)
let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr))
if area2 < 0.0
then
+#if DEBUG_LOG
+ printf "TWO area2=%f\n" area2
+#endif
area2 <- area2 + a2 * b2
area1 + area2
for i in 0..2 do
if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
then
- tanpts <- tanpts + 1;
- tanindex <- i;
-
+ tanpts <- tanpts + 1
+ tanindex <- i
+#if DEBUG_LOG
+ printf "tanindex=%d\n" tanindex
+#endif
+
if tanpts <> 1
then
-1.0
xint.[i] <- if xint.[i] < 0.0 then -a1 else a1
theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1)
+#if DEBUG_LOG
+ for k in 0..3 do
+ printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
+#endif
+
for j in 1..3 do
let tmp0 = theta.[j]
let tmp1 = xint.[j]
let tmp2 = yint.[j]
let mutable k = j - 1
+ let mutable k2 = 0
while k >= 0 do
if theta.[k] <= tmp0
then
+ k2 <- k + 1
k <- -1
else
theta.[k+1] <- theta.[k]
xint.[k+1] <- xint.[k]
yint.[k+1] <- yint.[k]
- k <- k - 1
+ k <- k - 1
+ k2 <- k + 1
+
+ theta.[k2] <- tmp0
+ xint.[k2] <- tmp1
+ yint.[k2] <- tmp2
- theta.[k+1] <- tmp0
- xint.[k+1] <- tmp1
- yint.[k+1] <- tmp2
+
+#if DEBUG_LOG
+ printf "AFTER sorting\n"
+ for k in 0..3 do
+ printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
+#endif
let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
if area5 < 0.0
then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
+#endif
area5 <- area5 + area_2
if area4 < 0.0
then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area4 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
+#endif
area4 <- area4 + area_2
if area3 < 0.0
then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area3 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
+#endif
area3 <- area3 + area_1
if area2 < 0.0
- then
+ then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
+#endif
area2 <- area2 + area_1
+#if DEBUG_LOG
+ printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5
+#endif
+
area1 + area2 + area3 + area4 + area5
quad ()
-open Types
-
-let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float =
- let { cx = h1; cy = k1; a = a1; b = b1; alpha = phi_1 } = e1
- let { cx = h2; cy = k2; a = a2; b = b2; alpha = phi_2 } = e2
+let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float =
+ 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
if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
then
let sinphi = sin phi_1
(h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
+#if DEBUG_LOG
+ printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
+#endif
+
let cosphi = cos phi_2r
let cosphi2 = cosphi ** 2.0
let sinphi = sin phi_2r
a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
|]
+#if DEBUG_LOG
+ for i in 0..4 do
+ printf "cy[%d]=%f\n" i cy.[i]
+#endif
+
let py = Array.zeroCreate<float> 5
let r = Array2D.zeroCreate<float> 3 5
for i in 0 .. 3 do
py.[4-i] <- cy.[i] / cy.[4]
py.[0] <- 1.0
+#if DEBUG_LOG
+ for i in 0..4 do
+ printf "py[%d]=%f\n" i py.[i]
+#endif
biquadroots py r
4
else
0
+#if DEBUG_LOG
+ 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
+#if DEBUG_LOG
+ printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
+#endif
|]
Array.sortInPlace ychk
+#if DEBUG_LOG
+ printf "nychk=%d\n" ychk.Length
+ for j in 0 .. ychk.Length - 1 do
+ printf "\t j=%d, ychk=%f\n" j ychk.[j]
+#endif
+
let nychk = Array.length ychk
let mutable nintpts = 0
let mutable i = 0
while returnValue = 0.0 && i < nychk do
+#if DEBUG_LOG
+ printf "------------->i=%d (nychk=%d)\n" i nychk
+#endif
+
if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
then
+#if DEBUG_LOG
+ printf "check intersecting points. nintps is %d" nintpts
+#endif
+
let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
let x2 = -x1
+#if DEBUG_LOG
+ printf "\tx1=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x1)= %f\n" x1 ychk.[i] a1 b1 (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff)
+ printf "\tx2=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x2) %f\n" x2 ychk.[i] a1 b1 (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff)
+#endif
+
if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
then
nintpts <- nintpts + 1
+#if DEBUG_LOG
+ printf "first if x1. acc nintps=%d\n" nintpts
+#endif
if nintpts > 4
then
returnValue <- -1.0
else
xint.[nintpts-1] <- x1
yint.[nintpts-1] <- ychk.[i]
+#if DEBUG_LOG
+ printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
+#endif
if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
then
nintpts <- nintpts + 1
+#if DEBUG_LOG
+ printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
+#endif
if nintpts > 4
then
returnValue <- -1.0
else
xint.[nintpts-1] <- x2
yint.[nintpts-1] <- ychk.[i]
+
+#if DEBUG_LOG
+ printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
+#endif
+
+#if DEBUG_LOG
+ else
+ printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
+#endif
i <- i + 1
match nintpts with
| 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
| 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
- | 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
- | INTERSECTION_POINT -> twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
+ | TANGENT_POINT ->
+#if DEBUG_LOG
+ printf "one point is tangent\n"
+#endif
+ nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
+
+ | INTERSECTION_POINT ->
+#if DEBUG_LOG
+ printf "check twointpts\n"
+#endif
+ twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| _ -> -1.0
-
projects / master-thesis.git / blobdiff