diff --git a/cue-ios/CueIOS/Views/EISView.swift b/cue-ios/CueIOS/Views/EISView.swift index 7cb8b15..e2e1234 100644 --- a/cue-ios/CueIOS/Views/EISView.swift +++ b/cue-ios/CueIOS/Views/EISView.swift @@ -285,46 +285,68 @@ struct EISView: View { with: .color(nyqColor)) } - // circle fit - if let fit = kasaCircleFit(points: points.map { (Double($0.zReal), Double(-$0.zImag)) }) { - let disc = fit.r * fit.r - fit.cy * fit.cy - if disc > 0 { - let sd = sqrt(disc) - let rs = fit.cx - sd - let rp = 2 * sd + // fit overlay + let fitColor = Color(red: 0.7, green: 0.3, blue: 0.9).opacity(0.5) + let fitPtColor = Color(red: 0.9, green: 0.4, blue: 0.9) - if rp > 0 { - let thetaR = atan2(-fit.cy, sd) - var thetaL = atan2(-fit.cy, -sd) - if thetaL < thetaR { thetaL += 2 * .pi } + if let result = fitNyquist(points: points.map { (Double($0.zReal), Double(-$0.zImag)) }) { + switch result { + case .circle(let fit): + let disc = fit.r * fit.r - fit.cy * fit.cy + if disc > 0 { + let sd = sqrt(disc) + let rs = fit.cx - sd + let rp = 2 * sd - let nArc = 120 - var arcPath = Path() - for i in 0...nArc { - let t = thetaR + (thetaL - thetaR) * Double(i) / Double(nArc) - let ax = fit.cx + fit.r * cos(t) - let ay = fit.cy + fit.r * sin(t) - let pt = CGPoint(x: lx(CGFloat(ax)), y: ly(CGFloat(ay))) - if i == 0 { arcPath.move(to: pt) } else { arcPath.addLine(to: pt) } + if rp > 0 { + let thetaR = atan2(-fit.cy, sd) + var thetaL = atan2(-fit.cy, -sd) + if thetaL < thetaR { thetaL += 2 * .pi } + + let nArc = 120 + var arcPath = Path() + for i in 0...nArc { + let t = thetaR + (thetaL - thetaR) * Double(i) / Double(nArc) + let ax = fit.cx + fit.r * cos(t) + let ay = fit.cy + fit.r * sin(t) + let pt = CGPoint(x: lx(CGFloat(ax)), y: ly(CGFloat(ay))) + if i == 0 { arcPath.move(to: pt) } else { arcPath.addLine(to: pt) } + } + context.stroke(arcPath, with: .color(fitColor), lineWidth: 1.5) + + let rsScr = CGPoint(x: lx(CGFloat(rs)), y: ly(0)) + let rpScr = CGPoint(x: lx(CGFloat(rs + rp)), y: ly(0)) + context.fill(Path(ellipseIn: CGRect(x: rsScr.x - 5, y: rsScr.y - 5, width: 10, height: 10)), + with: .color(fitPtColor)) + context.fill(Path(ellipseIn: CGRect(x: rpScr.x - 5, y: rpScr.y - 5, width: 10, height: 10)), + with: .color(fitPtColor)) + context.draw( + Text(String(format: "Rs=%.0f", rs)).font(.caption2).foregroundStyle(fitPtColor), + at: CGPoint(x: rsScr.x, y: rsScr.y + 14)) + context.draw( + Text(String(format: "Rp=%.0f", rp)).font(.caption2).foregroundStyle(fitPtColor), + at: CGPoint(x: rpScr.x, y: rpScr.y + 14)) } - context.stroke(arcPath, with: .color(Color(red: 0.7, green: 0.3, blue: 0.9).opacity(0.5)), - lineWidth: 1.5) - - // Rs and Rp markers - let fitPtColor = Color(red: 0.9, green: 0.4, blue: 0.9) - let rsScr = CGPoint(x: lx(CGFloat(rs)), y: ly(0)) - let rpScr = CGPoint(x: lx(CGFloat(rs + rp)), y: ly(0)) - context.fill(Path(ellipseIn: CGRect(x: rsScr.x - 5, y: rsScr.y - 5, width: 10, height: 10)), - with: .color(fitPtColor)) - context.fill(Path(ellipseIn: CGRect(x: rpScr.x - 5, y: rpScr.y - 5, width: 10, height: 10)), - with: .color(fitPtColor)) - context.draw( - Text(String(format: "Rs=%.0f", rs)).font(.caption2).foregroundStyle(fitPtColor), - at: CGPoint(x: rsScr.x, y: rsScr.y + 14)) - context.draw( - Text(String(format: "Rp=%.0f", rp)).font(.caption2).foregroundStyle(fitPtColor), - at: CGPoint(x: rpScr.x, y: rpScr.y + 14)) } + case .linear(let fit): + let xVals = points.map { CGFloat($0.zReal) }.filter { $0.isFinite } + guard let xMin = xVals.min(), let xMax = xVals.max() else { break } + let pad = (xMax - xMin) * 0.1 + let x0 = Double(xMin - pad) + let x1 = Double(xMax + pad) + let y0 = fit.slope * x0 + fit.yIntercept + let y1 = fit.slope * x1 + fit.yIntercept + var linePath = Path() + linePath.move(to: CGPoint(x: lx(CGFloat(x0)), y: ly(CGFloat(y0)))) + linePath.addLine(to: CGPoint(x: lx(CGFloat(x1)), y: ly(CGFloat(y1)))) + context.stroke(linePath, with: .color(fitColor), lineWidth: 1.5) + + let rsScr = CGPoint(x: lx(CGFloat(fit.rs)), y: ly(0)) + context.fill(Path(ellipseIn: CGRect(x: rsScr.x - 5, y: rsScr.y - 5, width: 10, height: 10)), + with: .color(fitPtColor)) + context.draw( + Text(String(format: "Rs=%.0f", fit.rs)).font(.caption2).foregroundStyle(fitPtColor), + at: CGPoint(x: rsScr.x, y: rsScr.y + 14)) } } } @@ -370,7 +392,7 @@ struct EISView: View { } } -// MARK: - Kasa circle fit (ported from plot.rs) +// MARK: - Nyquist fit struct CircleFitResult { let cx: Double @@ -378,6 +400,17 @@ struct CircleFitResult { let r: Double } +struct LinearFitResult { + let slope: Double + let yIntercept: Double + let rs: Double +} + +enum NyquistFitResult { + case circle(CircleFitResult) + case linear(LinearFitResult) +} + func kasaCircleFit(points: [(Double, Double)]) -> CircleFitResult? { let all = points.filter { $0.0.isFinite && $0.1.isFinite } guard all.count >= 4 else { return nil } @@ -463,6 +496,57 @@ private func kasaFitRaw(_ pts: [(Double, Double)]) -> (Double, Double, Double)? return (cx, cy, sqrt(rSq)) } +private func cumulativeTurning(_ pts: [(Double, Double)]) -> Double { + guard pts.count >= 3 else { return 0 } + var total = 0.0 + for i in 1..<(pts.count - 1) { + let (dx1, dy1) = (pts[i].0 - pts[i-1].0, pts[i].1 - pts[i-1].1) + let (dx2, dy2) = (pts[i+1].0 - pts[i].0, pts[i+1].1 - pts[i].1) + let cross = dx1 * dy2 - dy1 * dx2 + let dot = dx1 * dx2 + dy1 * dy2 + total += abs(atan2(cross, dot)) + } + return total +} + +private func fitLinear(_ pts: [(Double, Double)]) -> LinearFitResult? { + guard pts.count >= 2 else { return nil } + let n = Double(pts.count) + let sx = pts.map(\.0).reduce(0, +) + let sy = pts.map(\.1).reduce(0, +) + let sx2 = pts.map { $0.0 * $0.0 }.reduce(0, +) + let sxy = pts.map { $0.0 * $0.1 }.reduce(0, +) + let denom = n * sx2 - sx * sx + guard abs(denom) > 1e-20 else { return nil } + let slope = (n * sxy - sx * sy) / denom + let yInt = (sy - slope * sx) / n + let rs = abs(slope) > 1e-10 ? -yInt / slope : sx / n + return LinearFitResult(slope: slope, yIntercept: yInt, rs: rs) +} + +func fitNyquist(points: [(Double, Double)]) -> NyquistFitResult? { + let all = points.filter { $0.0.isFinite && $0.1.isFinite } + guard all.count >= 4 else { return nil } + + if cumulativeTurning(all) < 0.524 { + if let lin = fitLinear(all) { return .linear(lin) } + } + + if let circle = kasaCircleFit(points: points) { + let avgErr = all.map { p in + abs(sqrt((p.0 - circle.cx) * (p.0 - circle.cx) + + (p.1 - circle.cy) * (p.1 - circle.cy)) - circle.r) + }.reduce(0, +) / (Double(all.count) * circle.r) + if avgErr > 0.15 { + if let lin = fitLinear(all) { return .linear(lin) } + } + return .circle(circle) + } + + if let lin = fitLinear(all) { return .linear(lin) } + return nil +} + // MARK: - Canvas drawing helpers private func drawPolyline(context: GraphicsContext, points: [CGPoint], color: Color, width: CGFloat) {