sq = function(x) {
  return x * x;
}

nidx = _.map(d.conductors, String).indexOf(neutral)
pidx = _.map(d.conductors, String).indexOf(phase)

Yaddline = function(Y, Zseries, from, to) {
    var n = Zseries.x.length - 1;
    var Yseries = Zseries.inv();
    Y.setBlock([from,from], [from+n,from+n], Y.getBlock([from,from], [from+n,from+n]).add(Yseries));   // diagonal
    Y.setBlock([to,to],     [to+n,to+n],     Y.getBlock([to,to],     [to+n,to+n]).add(Yseries));
    Y.setBlock([from,to],   [from+n,to+n], Y.getBlock([from,to], [from+n,to+n]).sub(Yseries));   // off diagonal
    Y.setBlock([to,from],   [to+n,from+n], Y.getBlock([to,from], [to+n,from+n]).sub(Yseries));
    return Y;
}

Yaddshunt = function(Y, Zshunt, busnum) {
    var n = Zshunt.x.length;
    var Yshunt = Zshunt.inv();
    Y.setBlock([busnum,busnum], [busnum+n,busnum+n], Y.getBlock([busnum,busnum], [busnum+n,busnum+n]).add(Yshunt));   // diagonal
    return Y;
}

Yaddshort = function(Y, i, j) {
    var Ybig = 1e5
    Y.x[i][i] = Y.x[i][i] + Ybig
    Y.x[j][j] = Y.x[j][j] + Ybig
    Y.x[i][j] = Y.x[i][j] - Ybig
    Y.x[j][i] = Y.x[j][i] - Ybig
    return Y;
}

YaddlineRX = function(Y, R, X, i, j) {
    var k = sq(R) + sq(X)
    Y.x[i][i] = Y.x[i][i] + R/k
    Y.x[j][j] = Y.x[j][j] + R/k
    Y.x[i][j] = Y.x[i][j] - R/k
    Y.x[j][i] = Y.x[j][i] - R/k
    Y.y[i][i] = Y.y[i][i] - X/k
    Y.y[j][j] = Y.y[j][j] - X/k
    Y.y[i][j] = Y.y[i][j] + X/k
    Y.y[j][i] = Y.y[j][i] + X/k
    return Y;
}

numberOfSections = 50
subBus = 1
totalLength = totalLengthMi * 5.28  // ohms/kfeet
sectionLength = totalLength/numberOfSections  // kfeet
if (averageGround == "well") {
    rodlen = 200/3.28
    averageGround = rho/2/Math.PI/rodlen*(Math.log(4*rodlen/(8.6*0.0254) - 1))
} else if (isNaN(Number(averageGround))) {
    rodlen = 8/3.28
    averageGround = rho/2/Math.PI/rodlen*(Math.log(4*rodlen/(5/8*0.0254) - 1))
} else {
    averageGround = Number(averageGround)
}
Zgrnd = averageGround/(groundsPerKFeet * totalLength) * numberOfSections  // ohms (average ground per section)
d_ab = conductorSeparation  // feet (distance between the phase and neutral)

r_p = d.rac[pidx]/5.28  // ohms per 1000 feet
gmr_p = d.gmr[pidx]     // feet
r_n = d.rac[nidx]/5.28  // ohms per 1000 feet
gmr_n = d.gmr[nidx]     // feet

Z = numeric.t(numeric.identity(2), numeric.identity(2))
Z.x[0][0] = r_p + 0.01807
Z.y[0][0] = 0.0529 * math.log10(278.9 * math.sqrt(rho) / gmr_p)
Z.x[1][1] = r_n + 0.01807
Z.y[1][1] = 0.0529 * math.log10(278.9 * math.sqrt(rho) / gmr_n)
Z.x[0][1] = Z.x[1][0] = 0.01807
Z.y[0][1] = Z.y[1][0] = 0.0529 * math.log10(278.9 * math.sqrt(rho) / d_ab)
Zcond = Z.mul(sectionLength)   // actual ohms
numberOfConductors = 2

// preallocate the Ybus
n = (numberOfSections + 1) * numberOfConductors
Y = numeric.t(numeric.rep([n,n], 0.0), numeric.rep([n,n], 0.0))

// make the Ybus from the impedances
for (var i = 0; i < n - numberOfConductors - 1; i = i + numberOfConductors) {
    Y = Yaddline(Y, Zcond, i, i+numberOfConductors);
}

// add the shunt grounds
for (var i = 2*numberOfConductors - 1; i < n - 1; i = i + numberOfConductors) {
    Y.set([i, i], Y.get([i, i]).add(numeric.t(1/Zgrnd, 0.0)))
}

// add the sub ground
var i = 1
Y.set([i, i], Y.get([i, i]).add(numeric.t(1/Zsub, 0.0)))

// bond the phase wire to the neutral at fault location
faultNode = Math.floor(faultLoc * numberOfSections / 100) + 1
nfault = 2 * faultNode
Y = Yaddshort(Y, nfault - 2, nfault - 1)


// Add the source impedance
Y = YaddlineRX(Y, 0.0546 * systemVoltage/Math.sqrt(3) / faultI,  0.9985 * systemVoltage/Math.sqrt(3) / faultI, 0, 1)

// make current injections (Isrc)
Isrc = numeric.t(numeric.rep([n], 0.0), numeric.rep([n], 0.0))
Isrc.y[0] = -1000*faultI
Isrc.y[1] =  1000*faultI

// Find the voltages:
V = Y.inv().dot(Isrc)

// Find the voltage drops and conductor currents
Vdrops = numeric.t(numeric.rep([2,numberOfSections], 0.0), numeric.rep([2,numberOfSections], 0.0))
for (var i = 0; i < numberOfSections; i++) {
    Vdrops.set([0,i], V.get([i * 2    ]).sub(V.get([i * 2 + 2])))
    Vdrops.set([1,i], V.get([i * 2 + 1]).sub(V.get([i * 2 + 3])))
}
Ycond = Zcond.inv()
I = Ycond.dot(Vdrops)

vn = V.abs().x.filter(function(num,idx){ if( idx % 2 ) return num;})

println("NEV at the fault = " + math.format(vn[faultNode - 1]) + " V" + "  (" + math.format(vn[faultNode - 1] / systemVoltage / 1000 * Math.sqrt(3)) + " pu)")
println("Substation NEV = " + math.format(V.abs().x[1]) + " V")

    x = _.range(0, (numberOfSections+1) * sectionLength / 5.28, sectionLength / 5.28)
    seriesvn = _.zip(x,vn)
    $.plot($('#graph1'), [seriesvn]);
    
    seriesI = [{label: "Phase",
                data: _.zip(x,I.abs().x[0])},
               {label: "Neutral",
                data: _.zip(x,I.abs().x[1])},
               {label: "Ground",
                data: _.zip(x,I.getRow(0).add(I.getRow(1)).abs().x)}
               ]
    $.plot($('#graph2'), seriesI);