function renderNeutralInputs() {
    // Dynamically render neutral position inputs based on neutralCount
    const container = document.getElementById("neutralInputs");
    if (!container) return;
    const countEl = document.getElementById("neutralCount");
    const n = Math.max(1, parseInt(countEl && countEl.value ? countEl.value : 1, 10));
    let html = "";
    for (let i = 1; i <= n; i++) {
        const defaultX = i === 1 ? 10 : (i === 2 ? 40 : 10 + (i - 1) * 15);
        const defaultY = 71;
        html += `<label>Neutral ${i} position (x, y): </label>` +
                `<input type="text" id="n${i}x" name="n${i}x" value="${defaultX}" />, ` +
                `<input type="text" id="n${i}y" name="n${i}y" value="${defaultY}" /> (ft)<br />`;
    }
    container.innerHTML = html;
    setWarning("");
}

function setWarning(msg) {
    const el = document.getElementById("warning");
    if (el) {
        el.textContent = msg || "";
        el.style.display = msg ? "block" : "none";
    } else if (msg) {
        // Fallback alert if warning container is missing
        alert(msg);
    }
}

function calcImpedance() {
    const m_to_ft = (m) => m * 3.28084;
    const ft_to_m = (ft) => ft * 0.30479999;
    const mi_to_m = (mi) => mi * 1.60934 * 1000;
    const km_to_mi = (km) => km * 0.621371;
    const dist = (x1, y1, x2, y2) => Math.sqrt(Math.pow(x1 - x2, 2) + Math.pow(y1 - y2, 2));

    // System parameters
    const f = parseFloat(document.getElementById("frequency").value); // Hz

    // Conductor positions (feet)
    //const apos = [parseFloat(document.getElementById("ax".value)), parseFloat(document.getElementById("ay").value)];
    const apos = [parseFloat(document.getElementById("ax").value), parseFloat(document.getElementById("ay").value)];
    const bpos = [parseFloat(document.getElementById("bx").value), parseFloat(document.getElementById("by").value)];
    const cpos = [parseFloat(document.getElementById("cx").value), parseFloat(document.getElementById("cy").value)];
    // Neutral positions (feet) - arbitrary count
    const neutralCount = Math.max(1, parseInt(document.getElementById("neutralCount").value || "1", 10));
    const npos = [];
    for (let i = 1; i <= neutralCount; i++) {
        const xi = parseFloat(document.getElementById(`n${i}x`).value);
        const yi = parseFloat(document.getElementById(`n${i}y`).value);
        npos.push([xi, yi]);
    }
    // Validate no duplicate neutral positions (avoid log(0) singularities)
    const eps = 1e-9; // feet tolerance
    let dupMsg = "";
    for (let i = 0; i < neutralCount; i++) {
        for (let j = i + 1; j < neutralCount; j++) {
            const d = dist(npos[i][0], npos[i][1], npos[j][0], npos[j][1]);
            if (!isFinite(d) || isNaN(d)) { continue; }
            if (d < eps) {
                dupMsg = `Duplicate neutral positions detected: Neutral ${i+1} and Neutral ${j+1} share the same (x, y). Please adjust their positions.`;
                break;
            }
        }
        if (dupMsg) break;
    }
    if (dupMsg) {
        setWarning(dupMsg);
        return; // Abort calculation
    } else {
        setWarning("");
    }
    
    // Conductor properties
    const Ra = parseFloat(document.getElementById("Ra").value);           // Ohms per mile per conductor
    const Rb = Ra;
    const Rc = Ra;
    const GMRa = parseFloat(document.getElementById("GMRa").value);         // Geometric mean radius (feet)
    const GMRb = GMRa;
    const GMRc = GMRa;
    const ODa = parseFloat(document.getElementById("ODa").value);
    const bundle = parseFloat(document.getElementById("bundle").value);            // Conductors in bundle
    const bundle_spacing = parseFloat(document.getElementById("bundle_spacing").value);  // feet

    // Neutral conductor properties
    const Rn1 = parseFloat(document.getElementById("Rn1").value);            // Ohms per km per neutral conductor (assumed identical)
    const GMRn1 = parseFloat(document.getElementById("GMRn1").value);        // feet
    const ODn = parseFloat(document.getElementById("ODn").value);

    // values for dry earth
    let Dkkp = 1;
    if(document.getElementById("material").value == "seawater") {
        Dkkp = m_to_ft(38.25);  // feet
    } else if (document.getElementById("material").value == "swampyground") {
        Dkkp = m_to_ft(560);
    } else if (document.getElementById("material").value == "dampearth") {
        Dkkp = m_to_ft(850);
    } else if (document.getElementById("material").value == "dryearth") {
        Dkkp = m_to_ft(2690);
    } else if (document.getElementById("material").value == "pureslate") {
        Dkkp = m_to_ft(269000);
    } else if (document.getElementById("material").value == "sandstone") {
        Dkkp = m_to_ft(2690000);
    } else {
        Dkkp = m_to_ft(658.5*Math.pow((parseFloat(document.getElementById("resistivity").value)/f),0.5));
    }

    const Rkp = 9.869E-7 * f;
    const Rk = (Ra / mi_to_m(1)) / bundle; // convert to ohms per meter per phase
    const Rn = (Rn1 / 1000);                  // convert to ohms per meter per neutral

    let Dkkc = Math.pow(GMRa * Math.pow(bundle_spacing, (bundle - 1)), (1 / bundle));
    if (bundle >= 4){
        Dkkc = Dkkc * 1.091;
    }
    const Dkkg = GMRn1;

    const jwmu = math.multiply(math.complex(0,1), (2*Math.PI*f)*(2E-7));
    const Zaa = math.add(Rk + Rkp, math.multiply(jwmu, (math.log(Dkkp / Dkkc))));
    const Zbb = Zaa;
    const Zcc = Zaa;

    const Zab = math.add(Rkp, math.multiply(jwmu, (math.log(dist(apos[0], apos[1], bpos[0], bpos[1] - Dkkp) / dist(apos[0], apos[1], bpos[0], bpos[1])))));
    const Zbc = math.add(Rkp, math.multiply(jwmu, (math.log(dist(bpos[0], bpos[1], cpos[0], cpos[1] - Dkkp) / dist(bpos[0], bpos[1], cpos[0], cpos[1])))));
    const Zac = math.add(Rkp, math.multiply(jwmu, (math.log(dist(apos[0], apos[1], cpos[0], cpos[1] - Dkkp) / dist(apos[0], apos[1], cpos[0], cpos[1])))));

    // Build Zb (3 x N) and Zd (N x N)
    const Zb = [[], [], []];
    for (let k = 0; k < neutralCount; k++) {
        const nk = npos[k];
        const Za_n = math.add(Rkp, math.multiply(jwmu, (math.log(dist(apos[0], apos[1], nk[0], nk[1] - Dkkp) / dist(apos[0], apos[1], nk[0], nk[1])))));
        const Zb_n = math.add(Rkp, math.multiply(jwmu, (math.log(dist(bpos[0], bpos[1], nk[0], nk[1] - Dkkp) / dist(bpos[0], bpos[1], nk[0], nk[1])))));
        const Zc_n = math.add(Rkp, math.multiply(jwmu, (math.log(dist(cpos[0], cpos[1], nk[0], nk[1] - Dkkp) / dist(cpos[0], cpos[1], nk[0], nk[1])))));
        Zb[0][k] = Za_n;
        Zb[1][k] = Zb_n;
        Zb[2][k] = Zc_n;
    }

    const Zd = new Array(neutralCount).fill(null).map(()=> new Array(neutralCount).fill(null));
    for (let i = 0; i < neutralCount; i++) {
        for (let j = 0; j < neutralCount; j++) {
            if (i === j) {
                Zd[i][j] = math.add(Rn + Rkp, math.multiply(jwmu, (math.log(Dkkp / GMRn1))));
            } else {
                Zd[i][j] = math.add(Rkp, math.multiply(jwmu, (math.log(dist(npos[i][0], npos[i][1], npos[j][0], npos[j][1] - Dkkp) / dist(npos[i][0], npos[i][1], npos[j][0], npos[j][1])))));
            }
        }
    }

    const Zc = math.transpose(Zb);

    const Za = [
        [Zaa, Zab, Zac],
        [Zab, Zbb, Zbc],
        [Zac, Zbc, Zcc]
    ];

    const Zp = math.subtract(Za, math.multiply(math.multiply(Zb, math.inv(Zd)), Zc));
    
    const Zs = math.divide(math.add( math.add(Zp[0][0], Zp[1][1]), Zp[2][2]),3);
    const Zm = math.divide(math.add( math.add(Zp[0][1], Zp[0][2]), Zp[1][2]),3);

    console.log(Zs);
    console.log(Zm);
    const Zpsym = [
        [Zs, Zm, Zm],
        [Zm, Zs, Zm],
        [Zm, Zm, Zs]
    ];

    const a = math.complex(math.cos(120 * (Math.PI / 180)), math.sin(120 * (Math.PI / 180)));
    const Ainv = [
        [1, 1, 1],
        [1, a, math.pow(a, 2)],
        [1, math.pow(a, 2), a]
    ];
    const A = [
        [1, 1, 1],
        [1, math.pow(a, 2), a],
        [1, a, math.pow(a, 2)]
    ];

    const Zshat = math.multiply((1 / 3), math.multiply(math.multiply(Ainv, Zpsym), A));  // ohms per meter
    console.log(`Z0 = ${Zshat[0][0].re * 1000} + j${Zshat[0][0].im * 1000} ohm/km (zero sequence)`);
    console.log(`Z1 = ${Zshat[1][1].re * 1000} + j${Zshat[1][1].im * 1000} ohm/km (positive sequence)`);
    console.log(`Z2 = ${Zshat[2][2].re * 1000} + j${Zshat[2][2].im * 1000} ohm/km (negative sequence)`);

    document.getElementById("Z0").textContent = `${Zshat[0][0].re * 1000} + j${Zshat[0][0].im * 1000} ohm/km (zero sequence)`;
    document.getElementById("Z1").textContent = `${Zshat[1][1].re * 1000} + j${Zshat[1][1].im * 1000} ohm/km (positive sequence)`;
    document.getElementById("Z2").textContent = `${Zshat[2][2].re * 1000} + j${Zshat[2][2].im * 1000} ohm/km (negative sequence)`;


    const ra = ODa / (2 * 12); // inches to feet
    const rn = ODn / (2 * 12);

    let Dsc = Math.pow(ra * Math.pow(bundle_spacing, (bundle - 1)), (1 / bundle));
    if (bundle >= 4) {
        Dsc *= 1.091;
    }
    let Dscn = rn;

    const Paa = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((apos[1] * 2) / Dsc);
    const Pbb = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((bpos[1] * 2) / Dsc);
    const Pcc = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((cpos[1] * 2) / Dsc);

    const Pab = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((dist(apos[0], apos[1], bpos[0], (-1 * bpos[1]))) / (dist(apos[0], apos[1], bpos[0], bpos[1])));
    const Pac = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((dist(apos[0], apos[1], cpos[0], (-1 * cpos[1]))) / (dist(apos[0], apos[1], cpos[0], cpos[1])));
    const Pbc = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((dist(bpos[0], bpos[1], cpos[0], (-1 * cpos[1]))) / (dist(bpos[0], bpos[1], cpos[0], cpos[1])));

    // Build Pb (3 x N), Pc (N x 3), Pd (N x N)
    const Pb = [[], [], []];
    for (let k = 0; k < neutralCount; k++) {
        const nk = npos[k];
        const Pan = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((dist(apos[0], apos[1], nk[0], (-1 * nk[1]))) / (dist(apos[0], apos[1], nk[0], nk[1])));
        const Pbn = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((dist(bpos[0], bpos[1], nk[0], (-1 * nk[1]))) / (dist(bpos[0], bpos[1], nk[0], nk[1])));
        const Pcn = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((dist(cpos[0], cpos[1], nk[0], (-1 * nk[1]))) / (dist(cpos[0], cpos[1], nk[0], nk[1])));
        Pb[0][k] = Pan;
        Pb[1][k] = Pbn;
        Pb[2][k] = Pcn;
    }

    const Pd = new Array(neutralCount).fill(null).map(()=> new Array(neutralCount).fill(null));
    for (let i = 0; i < neutralCount; i++) {
        for (let j = 0; j < neutralCount; j++) {
            if (i === j) {
                Pd[i][j] = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((npos[i][1] * 2) / Dscn);
            } else {
                Pd[i][j] = (1 / (2 * Math.PI * 8.854E-12)) * Math.log((dist(npos[i][0], npos[i][1], npos[j][0], (-1 * npos[j][1]))) / (dist(npos[i][0], npos[i][1], npos[j][0], npos[j][1])));
            }
        }
    }

    const Pc = math.transpose(Pb);

    const Pa = [
        [Paa, Pab, Pac],
        [Pab, Pbb, Pbc],
        [Pac, Pbc, Pcc]
    ];
        
    const Cp = math.inv(math.subtract(Pa, math.multiply(math.multiply(Pb, math.inv(Pd)), Pc)));

    const Caa = math.divide(math.add( math.add(Cp[0][0], Cp[1][1]), Cp[2][2]),3);
    const Cab = math.divide(math.add( math.add(Cp[0][1], Cp[0][2]), Cp[1][2]),3);

    const Cpsym = [
        [Caa, Cab, Cab],
        [Cab, Caa, Cab],
        [Cab, Cab, Caa]
    ];

    //Cphat = np.multiply(1j*(2*np.pi*f)*1000,(1/3)*np.matmul(np.matmul(Ainv,Cpsym),A))
    // Cshat = (1/3)*np.matmul(np.matmul(Ainv,Cpsym),A)  //ohms per meter
    // 1j*2*np.pi*Cpsym
    const Cphat = math.multiply(math.multiply(math.complex(0,1), (2*Math.PI*f/3)), math.multiply(math.multiply(Ainv, Cpsym), A));
    console.log(Cphat);
    console.log(`Y0 = ${Cphat[0][0].im * 1000} ohm/km (zero sequence)`);
    console.log(`Y1 = ${Cphat[1][1].im * 1000} ohm/km (positive sequence)`);
    console.log(`Y2 = ${Cphat[2][2].im * 1000} ohm/km (negative sequence)`);

    document.getElementById("Y0").textContent = `j${Cphat[0][0].im * 1000} ohm/km (zero sequence)`;
    document.getElementById("Y1").textContent = `j${Cphat[1][1].im * 1000} ohm/km (positive sequence)`;
    document.getElementById("Y2").textContent = `j${Cphat[2][2].im * 1000} ohm/km (negative sequence)`;

}

function dropdownChangeHandler() {
    if(document.getElementById("material").value=="custom") {
        document.getElementById("resistivitylabel").style="display:inline;";
        document.getElementById("resistivity").style="display:inline;";
        document.getElementById("resistivityunits").style="display:inline;";
    } else {
        document.getElementById("resistivitylabel").style="display:none;";
        document.getElementById("resistivity").style="display:none;";
        document.getElementById("resistivityunits").style="display:none;";
    }
}

// Initialize dynamic fields on load
if (document.readyState === 'loading') {
    document.addEventListener('DOMContentLoaded', renderNeutralInputs);
} else {
    renderNeutralInputs();
}