import type { PileCapGeometry, MaterialProperties, LoadCombinationResult, ReinforcementConfig } from '../types/pileCap';

export interface ColumnDesignResult {
    capacityX: number; // k-ft
    capacityY: number; // k-ft
    demandX: number; // k-ft
    demandY: number; // k-ft
    axialLoad: number; // kips
    ratio: number;
    pass: boolean;
    interactionPointsX: { m: number, p: number }[];
    interactionPointsY: { m: number, p: number }[];
}

export function performColumnDesign(
    geometry: PileCapGeometry,
    materials: MaterialProperties,
    reinforcement: ReinforcementConfig,
    load: LoadCombinationResult
): ColumnDesignResult {
    const { columnLength, columnWidth } = geometry;
    const { columnBarSize, columnBarCount } = reinforcement;

    // Parse bar size
    const barSizeNum = parseInt(columnBarSize.replace('#', '')) || 9;
    const barDiameter = barSizeNum / 8;
    const barArea = Math.PI * Math.pow(barDiameter / 2, 2);
    const totalSteelArea = columnBarCount * barArea;

    // Generate Interaction Diagram for X-Axis (Bending about Y)
    // Column dimension h = columnLength (X-dir)
    // Width b = columnWidth (Z-dir)
    // Bars distributed: Assume equal distribution on all 4 sides for simplicity
    // Or simplified: Smear steel as a tube? Or discrete bars?
    // Let's model discrete bars.

    // Distribute bars
    // If count = 8, maybe 3 on each face (corners shared)? 
    // Top/Bottom rows: N_x bars
    // Side columns: N_y bars
    // Total = 2*Nx + 2*Ny - 4 = Count
    // Assume square arrangement if possible, or proportional to dimensions.

    // For simplicity in this checkpoint, let's assume bars are placed in a single layer around the perimeter.
    // We will generate points for the P-M interaction curve.

    const pointsX = generateInteractionDiagram(
        columnLength,
        columnWidth,
        totalSteelArea,
        materials.fc,
        materials.fy,
        2.5 // cover to center of bars
    );

    const pointsY = generateInteractionDiagram(
        columnWidth,
        columnLength,
        totalSteelArea,
        materials.fc,
        materials.fy,
        2.5
    );

    // Find Capacity for given Axial Load
    // Interpolate M capacity from P
    const Pu = load.axialLoad;
    const Mux = Math.abs(load.momentX); // k-ft
    const Muy = Math.abs(load.momentY); // k-ft

    const phiMnX = interpolateMomentCapacity(pointsX, Pu);
    const phiMnY = interpolateMomentCapacity(pointsY, Pu);

    // Biaxial Check (Bresler Load Contour Method)
    // (Mux / phiMnx)^alpha + (Muy / phiMny)^alpha <= 1.0
    // alpha depends on bar arrangement, usually 1.15 to 1.5. 
    // Conservative: alpha = 1.0 (linear)
    // ACI PCA method suggests alpha based on P/Po.
    // Let's use alpha = 1.0 for simplicity and conservatism.

    let ratio = 0;
    if (phiMnX > 0) ratio += Mux / phiMnX;
    if (phiMnY > 0) ratio += Muy / phiMnY;

    // If axial load exceeds max compression, ratio is infinite
    const maxP = pointsX[0].p; // Pure compression
    if (Pu > maxP) ratio = 999;

    return {
        capacityX: phiMnX,
        capacityY: phiMnY,
        demandX: Mux,
        demandY: Muy,
        axialLoad: Pu,
        ratio,
        pass: ratio <= 1.0,
        interactionPointsX: pointsX,
        interactionPointsY: pointsY
    };
}

function generateInteractionDiagram(
    h: number, // Dimension perpendicular to axis of bending (depth)
    b: number, // Dimension parallel to axis of bending (width)
    As: number,
    fc: number,
    fy: number,
    d_prime: number // distance from face to centroid of steel
): { m: number, p: number }[] {
    // Detailed Interaction Diagram Generation using Strain Compatibility
    // Iterate neutral axis depth 'c' from slightly larger than h (pure compression) down to 0 (pure tension)

    const points: { m: number, p: number }[] = [];
    const d = h - d_prime; // Effective depth
    const Es = 29000; // ksi
    const beta1 = fc <= 4000 ? 0.85 : Math.max(0.65, 0.85 - 0.05 * (fc - 4000) / 1000);

    // Steel layers (Simplified: 2 layers, top and bottom for this direction)
    // For a more accurate diagram, we should model the actual bar distribution.
    // Let's assume:
    // - As/2 at d (Tension/Compression)
    // - As/2 at d_prime (Compression/Tension)
    // This is a "two-face" simplification. For 4-face distributed bars, it's more complex.
    // Given the user's request for a "column graph", a 2-face model is often a reasonable approximation for checking,
    // but for a real design tool, we should ideally distribute.
    // Let's try to distribute: 3 layers? Top, Middle, Bottom.
    // If we assume uniform distribution around perimeter:
    // Top layer: As/4
    // Bottom layer: As/4
    // Side layers: As/2 distributed?
    // Let's stick to the 2-layer simplification for robustness in this step, or a 3-layer if easy.
    // 2-layer is standard for "approximate" curves unless we have exact bar coordinates.
    // Let's use 2 layers (Top and Bottom) which is conservative for bending if side bars are ignored for moment but included for axial.
    // Actually, including side bars increases moment capacity.
    // Let's stick to: Total As distributed as: 50% at top, 50% at bottom. This is a common simplification.

    const As_layer = As / 2;

    // 1. Pure Compression (Point 0)
    // P0 = 0.85 * fc * (Ag - As) + As * fy
    const Ag = h * b;
    const P0 = (0.85 * fc * (Ag - As) + As * fy) / 1000; // kips
    const phi_c_tied = 0.65;
    const alpha = 0.80;
    const phiPn_max = phi_c_tied * alpha * P0;

    // Add max axial point (with moment = 0)
    points.push({ m: 0, p: phiPn_max });

    // Iterate c from 1.5h down to 0.1h
    const steps = 20;
    const c_start = h * 1.1;
    const c_end = h * 0.1;
    const c_step = (c_start - c_end) / steps;

    for (let c = c_start; c >= c_end; c -= c_step) {
        // Strain diagram
        // eps_c = 0.003 at extreme compression fiber
        // eps_s = eps_c * (d - c) / c

        // Concrete compression block
        let a = beta1 * c;
        if (a > h) a = h; // Cap at h
        const Cc = 0.85 * fc * a * b; // lbs

        // Steel forces
        // Layer 1 (Compression face): at d_prime
        const eps_s1 = 0.003 * (c - d_prime) / c;
        const fs1 = Math.min(Math.max(eps_s1 * Es, -fy), fy); // ksi
        const Fs1 = As_layer * fs1 * 1000; // lbs (Compression is positive here for force sum?)
        // Wait, let's use consistent sign convention. Compression positive.

        // Layer 2 (Tension face): at d
        const eps_s2 = 0.003 * (c - d) / c; // Negative if tension
        // eps_s2 is strain at tension steel.
        // If c < d, (c-d) is negative -> tension.
        const fs2 = Math.min(Math.max(eps_s2 * Es, -fy), fy);
        const Fs2 = As_layer * fs2 * 1000; // lbs

        // Nominal Axial Strength Pn
        const Pn = (Cc + Fs1 + Fs2) / 1000; // kips

        // Nominal Moment Strength Mn (about plastic centroid = h/2)
        const M_Cc = Cc * (h / 2 - a / 2);
        const M_Fs1 = Fs1 * (h / 2 - d_prime);
        const M_Fs2 = Fs2 * (h / 2 - d); // d is > h/2, so this term handles sign correctly if Fs2 is signed

        const Mn = (M_Cc + M_Fs1 + M_Fs2) / 12000; // k-ft

        // Determine phi
        // Based on Net Tensile Strain (eps_t) in extreme tension steel
        // Here eps_t = abs(eps_s2) if eps_s2 is tension (negative)
        let phi = 0.65;
        const eps_t = -eps_s2; // Tension is negative in my calc above, so negate to get magnitude of tension strain

        if (eps_t <= 0.002) {
            phi = 0.65;
        } else if (eps_t >= 0.005) {
            phi = 0.90;
        } else {
            phi = 0.65 + 0.25 * (eps_t - 0.002) / 0.003;
        }

        // Cap Pn at Pn_max? No, the curve goes up to P0, but we slice it at Pn_max usually.
        // But for the curve generation, we plot the theoretical PhiPn, then cap it visually or logic-wise.
        // The logic `interpolateMomentCapacity` handles the check against max P.
        // So we just store the point.

        // Only add if Pn is positive (compression) or slightly tension?
        // Interaction diagram usually shows compression zone.
        if (Pn > -P0 * 0.1) { // Allow some tension for completeness
            points.push({ m: Mn * phi, p: Pn * phi });
        }
    }

    // Pure Tension
    // Pnt = As * fy
    const Pnt = -As * fy / 1000; // kips (Tension)
    const phi_t = 0.9;
    points.push({ m: 0, p: Pnt * phi_t });

    // Filter points to ensure we don't exceed Pn_max (horizontal cut)
    // Actually, we should keep the theoretical curve and let the renderer or checker handle the cap.
    // But standard diagrams often show the flat top.
    // Let's clamp the first few points if they exceed phiPn_max

    return points.map(pt => ({
        m: pt.m,
        p: Math.min(pt.p, phiPn_max)
    })).sort((a, b) => b.p - a.p); // Sort by P descending
}

function interpolateMomentCapacity(points: { m: number, p: number }[], Pu: number): number {
    // Find points bounding Pu
    // Points are sorted by P descending (Max P -> 0)

    if (Pu > points[0].p) return 0; // Exceeds max axial
    if (Pu < 0) return points[points.length - 1].m; // Tension? Assume Pure Moment for small tension for now

    for (let i = 0; i < points.length - 1; i++) {
        const p1 = points[i];
        const p2 = points[i + 1];

        if (Pu <= p1.p && Pu >= p2.p) {
            // Interpolate M
            const ratio = (Pu - p2.p) / (p1.p - p2.p);
            return p2.m + ratio * (p1.m - p2.m);
        }
    }

    return 0;
}