// Translation from https://github.com/chraibi/EEOver.
module ParasitemiaCore.EEOver

open System

let private EPS = 1.0e-7

let inline private ellipse2tr (x : float) (y : float) (aa : float) (bb : float) (cc : float) (dd : float) (ee : float) (ff : float) : float =
    aa * x * x + bb * x * y + cc * y * y + dd * x + ee * y + ff

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) =
    let a1b1 = a1 * b1
    let a2b2 = a2 * b2
    let area_1 = Math.PI * a1b1
    let area_2 = Math.PI * a2b2
    let relsize = a1b1 - a2b2

    if relsize > 0.0 then
        if (h2_tr * h2_tr) / (a1 * a1) + (k2_tr * k2_tr) / (b1 * b1) < 1.0 then
            area_2
        else
            0.0

    elif relsize < 0.0 then
        if ff < 0.0 then
            area_1
        else
            0.0

    else
        if abs (h1 - h2) < EPS && abs (k1 - k2) < EPS && abs (area_1 - area_2) < EPS then
            area_1
        else
            0.0

type private PointType = TANGENT_POINT | INTERSECTION_POINT

let private istanpt (x : float) (y : float) (a1 : float) (b1 : float) (aa : float) (bb : float) (cc : float) (dd : float) (ee : float) (ff : float) : PointType =
    let x =
        if abs x > a1 then
            if x < 0.0 then -a1 else a1
        else
            x

    let theta =
        if y < 0.0 then
            2.0 * Math.PI - acos (x / a1)
        else
            acos (x / a1)

    let eps_radian = 0.1

    let x1 = a1 * cos (theta + eps_radian)
    let y1 = b1 * sin (theta + eps_radian)
    let x2 = a1 * cos (theta - eps_radian)
    let y2 = b1 * sin (theta - eps_radian)

    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

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) =
    if abs x.[0] > a1 then
        x.[0] <- if x.[0] < 0.0 then -a1 else a1

    let mutable theta1 =
        if y.[0] < 0.0 then
            2.0 * Math.PI - acos (x.[0] / a1)
        else
            acos (x.[0] / a1)

    if abs x.[1] > a1 then
        x.[1] <- if x.[1] < 0.0 then -a1 else a1

    let mutable theta2 =
        if y.[1] < 0.0 then
            2.0 * Math.PI - acos (x.[1] / a1)
        else
            acos (x.[1] / a1)

    if theta1 > theta2 then
        let tmp = theta1
        theta1 <- theta2
        theta2 <- tmp

    let xmid = a1 * cos ((theta1 + theta2) / 2.0)
    let ymid = b1 * sin ((theta1 + theta2) / 2.0)

    if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0 then
        let tmp = theta1
        theta1 <- theta2
        theta2 <- tmp

    if theta1 > theta2 then
        theta1 <- theta1 - 2.0 * Math.PI

    let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0

    let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0]))

    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 sinphi = sin (phi_1 - phi_2)

    let mutable x1_tr = (x.[0] - h2_tr) * cosphi + (y.[0] - k2_tr) * -sinphi
    let mutable y1_tr = (x.[0] - h2_tr) * sinphi + (y.[0] - k2_tr) * cosphi
    let mutable x2_tr = (x.[1] - h2_tr) * cosphi + (y.[1] - k2_tr) * -sinphi
    let mutable y2_tr = (x.[1] - h2_tr) * sinphi + (y.[1] - k2_tr) * cosphi

    if abs x1_tr > a2 then
        x1_tr <- if x1_tr < 0.0 then -a2 else a2

    if y1_tr < 0.0 then
        theta1 <- 2.0 * Math.PI - acos (x1_tr / a2)
    else
        theta1 <- acos (x1_tr / a2)

    if abs x2_tr > a2 then
        x2_tr <- if x2_tr < 0.0 then -a2 else a2

    if y2_tr < 0.0 then
        theta2 <- 2.0 * Math.PI - acos (x2_tr / a2)
    else
        theta2 <- acos (x2_tr / a2)

    if theta1 > theta2 then
        let tmp = theta1
        theta1 <- theta2
        theta2 <- tmp

    let xmid = a2 * cos ((theta1 + theta2) / 2.0)
    let ymid = b2 * sin ((theta1 + theta2) / 2.0)

    let cosphi = cos (phi_2 - phi_1)
    let sinphi = sin (phi_2 - phi_1)
    let xmid_rt = xmid * cosphi + ymid * -sinphi + h2_tr
    let ymid_rt = xmid * sinphi + ymid * cosphi + k2_tr

    if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0 then
        let tmp = theta1
        theta1 <- theta2
        theta2 <- tmp

    if theta1 > theta2 then
        theta1 <- theta1 - 2.0 * Math.PI

    let trsign = if theta2 - theta1 > Math.PI then 1.0 else -1.0

    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

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 =
    let mutable tanpts = 0
    let mutable tanindex = 0
    for i = 0 to 2 do
        if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT then
            tanpts <- tanpts + 1
            tanindex <- i
#if DEBUG_LOG
    printf "tanindex=%d\n" tanindex
#endif

    if tanpts <> 1 then
        -1.0
    else
        match tanindex with
        | 0 ->
            xint.[0] <- xint.[2]
            yint.[0] <- yint.[2]
        | 1 ->
            xint.[1] <- xint.[2]
            yint.[1] <- yint.[2]
        | _ ->
            ()
        twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff

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 =
    let a1b1 = a1 * b1
    let a2b2 = a2 * b2
    let area_1 = Math.PI * a1b1
    let area_2 = Math.PI * a2b2

    let theta = Array.zeroCreate 4

    for i = 0 to 3 do
        if abs xint.[i] > a1 then
            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 = 0 to 3 do
        printf "k=%d:  Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
#endif

    for j = 1 to 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
                k2 <- k + 1

        theta.[k2] <- tmp0
        xint.[k2] <- tmp1
        yint.[k2] <- tmp2


#if DEBUG_LOG
    printf "AFTER sorting\n"
    for k = 0 to 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]))

    let cosphi = cos (phi_1 - phi_2)
    let sinphi = sin (phi_1 - phi_2)

    let theta_tr = Array.zeroCreate 4
    let xint_tr = Array.zeroCreate 4
    let yint_tr = Array.zeroCreate 4

    for i = 0 to 3 do
        xint_tr.[i] <- (xint.[i] - h2_tr) * cosphi + (yint.[i] - k2_tr) * -sinphi
        yint_tr.[i] <- (xint.[i] - h2_tr) * sinphi + (yint.[i] - k2_tr) * cosphi

        if abs xint_tr.[i] > a2 then
            xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2

        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)

    let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0)
    let ymid = b1 * sin ((theta.[0] + theta.[1]) / 2.0)

    let mutable area2, area3, area4, area5 = 0.0, 0.0, 0.0, 0.0

    if ellipse2tr xmid ymid aa bb cc dd ee ff < 0.0 then
        area2 <- 0.5 * (a1b1 * (theta.[1] - theta.[0]) - abs (xint.[0] * yint.[1] - xint.[1] * yint.[0]))
        area3 <- 0.5 * (a1b1 * (theta.[3] - theta.[2]) - abs (xint.[2] * yint.[3] - xint.[3] * yint.[2]))
        area4 <- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs (xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1]))

        if theta_tr.[3] > theta_tr.[0] then
            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]))
        else
            area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
    else
        area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs (xint.[1] * yint.[2] - xint.[2] * yint.[1]))
        area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3]))
        area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
        area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))

    if area5 < 0.0 then
#if DEBUG_LOG
        printf "\n\t\t-------------> area5 is negative (%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 negative (%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 negative (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
#endif
        area3 <- area3 + area_1

    if area2 < 0.0 then
#if DEBUG_LOG
        printf "\n\t\t-------------> area2 is negative (%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

let private quadroots (p : float[]) (r : float[,]) =
    let mutable b = -p.[1] / (2.0 * p.[0])
    let c = p.[2] / p.[0]
    let mutable d = b * b - c

    if d >= 0.0 then
        if b > 0.0 then
            b <- sqrt d + b
            r.[1, 2] <- b
        else
            b <- -sqrt d + b
            r.[1, 2] <- b
        r.[1, 1] <- c / b
        r.[2, 1] <- 0.0
        r.[2, 2] <- 0.0
    else
        d <- sqrt -d
        r.[2, 1] <- d
        r.[2, 2] <- -d
        r.[1, 1] <- b
        r.[1, 2] <- b

let private cubicroots (p : float[]) (r : float[,]) =
    if p.[0] <> 1.0 then
        for k = 1 to 3 do
            p.[k] <- p.[k] / p.[0]
        p.[0] <- 1.0
    let s = p.[1] / 3.0
    let mutable t = s * p.[1]
    let mutable b = 0.5 * (s * (t / 1.5 - p.[2]) + p.[3])
    t <- (t - p.[2]) / 3.0
    let mutable c = t * t * t
    let mutable d = b * b - c

    if d >= 0.0 then
        d <- (sqrt d + abs b) ** (1.0 / 3.0)
        if d <> 0.0 then
            if b > 0.0 then
                b <- -d
            else
                b <- d
            c <- t / b
        d <- sqrt 0.75 * (b - c)
        r.[2, 2] <- d
        b <- b + c
        c <- -0.5 * b - s
        r.[1, 2] <- c
        if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0 then
            r.[1, 1] <- c
            r.[2, 1] <- -d
            r.[1, 3] <- b - s
            r.[2, 3] <- 0.0
        else
            r.[1, 1] <- b - s
            r.[2, 1] <- 0.0
            r.[1, 3] <- c
            r.[2, 3] <- -d
    else
        if b = 0.0 then
            d <- (atan 1.0) / 1.5
        else
            d <- atan ((sqrt -d) / (abs b)) / 3.0

        if b < 0.0 then
            b <- 2.0 * (sqrt t)
        else
            b <- -2.0 * (sqrt t)

        c <- (cos d) * b
        t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c
        d <- -t - c - s
        c <- c - s
        t <- t - s

        if abs c > abs t then
            r.[1, 3] <- c
        else
            r.[1, 3] <- t
            t <- c

        if abs d > abs t then
            r.[1, 2] <- d
        else
            r.[1, 2] <- t
            t <- d

        r.[1, 1] <- t
        for k = 1 to 3 do
            r.[2, k] <- 0.0

let inline private biquadroots (p : float[]) (r : float[,]) =
    if p.[0] <> 1.0 then
        for k = 1 to 4 do
            p.[k] <- p.[k] / p.[0]
        p.[0] <- 1.0
    let e = 0.25 * p.[1]
    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])
    a <- a - d

    let mutable quadExecuted = false
    let inline quad () =
        if not quadExecuted then
            p.[2] <- c / b
            quadroots p r
            for k = 1 to 2 do
                for j = 1 to 2 do
                    r.[j, k+2] <- r.[j, k]
            p.[1] <- -p.[1]
            p.[2] <- b
            quadroots p r
            for k = 1 to 4 do
                r.[1,k] <- r.[1,k] - e
            quadExecuted <- true

    p.[1] <- 0.5 * a
    p.[2] <- (p.[1] * p.[1] - c) * 0.25
    p.[3] <- b * b / -64.0
    if p.[3] < 0.0 then
        cubicroots p r
        let mutable k = 1
        while k < 4 do
            if r.[2, k] = 0.0 && r.[1, k] > 0.0 then
                d <- r.[1, k] * 4.0
                a <- a + d
                if a >= 0.0 && b >= 0.0 then
                    p.[1] <- sqrt d
                elif a <= 0.0 && b <= 0.0 then
                    p.[1] <- sqrt d
                else
                    p.[1] <- -(sqrt d)
                b <- 0.5 * (a + b / p.[1])
                quad ()
                k <- 4
            k <- k + 1

    if not quadExecuted && p.[2] < 0.0 then
        b <- sqrt c
        d <- b + b - a
        p.[1] <- 0.0
        if d > 0.0 then
            p.[1] <- sqrt d
    elif not quadExecuted then
        if p.[1] > 0.0 then
            b <- (sqrt p.[2]) * 2.0 + p.[1]
        else
            b <- -(sqrt p.[2]) * 2.0 + p.[1]

        if b <> 0.0 then
            p.[1] <- 0.0
        else
            for k = 1 to 4 do
                r.[1, k] <- -e
                r.[2, k] <- 0.0
            quadExecuted <- true

    quad ()

/// <summary>
/// Return a tuple (area, x intersections, y intersections).
/// </summary>
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
        None
    else
        let phi_1 = phi_1 % Math.PI
        let phi_2 = phi_2 % Math.PI
        let h2_tr, k2_tr, phi_2r =
            let cosphi = cos phi_1
            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
        let sinphi2 = sinphi ** 2.0
        let cosphisinphi = 2.0 * cosphi * sinphi
        let a22 = a2 ** 2.0
        let b22 = b2 ** 2.0
        let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22
        let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22
        let tmp2 = cosphi * h2_tr + sinphi * k2_tr
        let tmp3 = sinphi * h2_tr - cosphi * k2_tr

        let aa = cosphi2 / a22 + sinphi2 / b22
        let bb = cosphisinphi / a22 - cosphisinphi / b22
        let cc = sinphi2 / a22 + cosphi2 / b22
        let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1
        let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1
        let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0

        let cy = [|
            (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff)
            2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff)
            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)
            2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee))
            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 = 0 to 4 do
            printf "cy[%d]=%f\n" i cy.[i]
#endif

        let py = Array.zeroCreate<float> 5
        let r = Array2D.zeroCreate<float> 3 5

        let nroots =
            if abs cy.[4] > EPS then
                for i = 0 to 3 do
                    py.[4-i] <- cy.[i] / cy.[4]
                py.[0] <- 1.0
#if DEBUG_LOG
                for i to 0 to 4 do
                    printf "py[%d]=%f\n" i py.[i]
#endif
                biquadroots py r
                4

            elif abs cy.[3] > EPS then
                for i = 0 to 2 do
                    py.[3 - i] <- cy.[i] / cy.[3]
                py.[0] <- 1.0
                cubicroots py r
                3

            elif abs cy.[2] > EPS then
                for i = 0 to 1 do
                    py.[2-i] <- cy.[i] / cy.[2]
                py.[0] <- 1.0
                quadroots py r
                2

            elif abs cy.[1] > EPS then
                r.[1, 1] <- -cy.[0] / cy.[1]
                r.[2, 1] <- 0.0
                1

            else
                0

#if DEBUG_LOG
        printf "nroots = %d\n" nroots
#endif

        let ychk = Array.init nroots (fun _ -> Double.MaxValue)
        let mutable nychk = 0
        for i = 1 to 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
#endif
        Array.sortInPlace ychk

#if DEBUG_LOG
        printf "nychk=%d\n" ychk.Length
        for j = 0 to ychk.Length - 1 do
            printf "\t j=%d, ychk=%f\n" j ychk.[j]
#endif

        let mutable nintpts = 0

        let xint = Array.zeroCreate 4
        let yint = Array.zeroCreate 4

        let mutable returnValue = 0.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 "\n\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

        if returnValue = -1.0 then
            None
        else
            let area =
                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 ->
#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

            if area = -1.0 then
                None
            elif nintpts = 0 then
                Some (float32 area, [||], [||])
            else
                let xTransform : float32[] = Array.zeroCreate nintpts
                let yTransform : float32[] = Array.zeroCreate nintpts
                for i = 0 to (nintpts - 1) do
                    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)