if (harmonicSelect == "Fundamental") {
    var harmonic = 1
} else {
    var harmonic = 3
}


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

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

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;
}

numberOfSections = 50
subBus = 1
totalLength = totalLengthMi * 5.28  // ohms/kfeet
sectionLength = totalLength/numberOfSections  // kfeet
Zgrnd = averageGround/(groundsPerKFeet * totalLength) * numberOfSections  // ohms (average ground per section)
d_ab = 6  // feet (distance between the phase and neutral)

r = 0.0001  // ohms per 1000 feet
gmr = 0.04  // 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 + 0.01807 * harmonic
Z.y[0][0] = harmonic * 0.0529 * math.log10(278.9 * math.sqrt(rho/harmonic) / gmr)
Z.x[1][1] = r_n + 0.01807 * harmonic
Z.y[1][1] = harmonic * 0.0529 * math.log10(278.9 * math.sqrt(rho/harmonic) / gmr_n)
Z.x[0][1] = Z.x[1][0] = 0.01807 * harmonic
Z.y[0][1] = Z.y[1][0] = harmonic * 0.0529 * math.log10(278.9 * math.sqrt(rho/harmonic) / 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)))

// ground the phase wire at the sub
var i = 0
var j = 1
Ylarge = numeric.t(100000,0.0)
Y.set([i, i], Y.get([i, i]).add(Ylarge))
Y.set([j, j], Y.get([j, j]).add(Ylarge))
Y.set([i, j], Y.get([i, j]).sub(Ylarge))
Y.set([j, i], Y.get([j, i]).sub(Ylarge))

// make current injections (Isrc)
Isrc = numeric.t(numeric.rep([n], 0.0), numeric.rep([n], 0.0))
// inject evenly distributed loads
if (loadType == "Evenly distributed") {
    for (var i = numberOfConductors*subBus; i < n; i = i + numberOfConductors) {
        // phase A
        Isrc.set([i],   numeric.t( UnbalanceI/(numberOfSections-subBus), 0.0)) // amperes
        // neutral
        Isrc.set([i+1], numeric.t(-UnbalanceI/(numberOfSections-subBus), 0.0)) // amperes
    }
} else if (loadType == "Lumped at end") {
    // inject isolated load at the end
    Isrc.set([n-2], numeric.t(-UnbalanceI, 0.0)) // amperes
    Isrc.set([n-1], numeric.t( UnbalanceI, 0.0)) // amperes
} else {  // substation injection
    Isrc.set([1], numeric.t(UnbalanceI, 0.0)) // amperes
}

// 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("Substation NEV = " + math.format(V.abs().x[1]) + " V")
println("Maximum NEV = " + math.format(_.max(vn)) + " 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);