Y.addline = function(Y, Zseries, from, to) {
  n = NROW(Zseries)
  fromNodes = from + 1:n - 1
  toNodes = to + 1:n - 1
  Yseries = solve(Zseries)
  Y[fromNodes,fromNodes] = Y[fromNodes,fromNodes] + Yseries # diagonal
  Y[toNodes,toNodes] = Y[toNodes,toNodes] + Yseries
  Y[fromNodes,toNodes] = Y[fromNodes,toNodes] - Yseries # off diagonal
  Y[toNodes,fromNodes] = Y[toNodes,fromNodes] - Yseries
  Y
}

Y.addshunt = function(Y, Zshunt, busnum) {
  n = NROW(Zshunt)
  fromNodes = busnum:(busnum+n-1)
  Yshunt = solve(Zshunt)
  Y[fromNodes,fromNodes] = Y[fromNodes,fromNodes] + Yshunt # diagonal
  Y
}
load(RpadBaseFile("AA_conductors.RData"))
HTMLon()
cat("Neutral conductors: ")
HTMLselect("neutralConductor1", a$label, default=16)
HTMLselect("neutralConductor2", c("NONE",a$label))

numberOfSections = 50
subBus = 20
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)

aa.n1 = a[a$label==neutralConductor1,]

r = 0.0001 # ohms per 1000 feet
gmr = 0.04 # feet
r.n = aa.n1$rac/5.28 # ohms per 1000
gmr.n = aa.n1$gmr

if (neutralConductor2 != "NONE") {
  aa.n2 = a[a$label==neutralConductor2,]
  r.n2 = aa.n2$rac/5.28 # ohms per 1000
  gmr.n2 = aa.n2$gmr
  Zcond = array(0i,c(3,3)) # ohms per 1000 feet
  r.n2 = aa.n2$rac/5.28 # ohms per 1000
  gmr.n2 = aa.n2$gmr
  Zcond[3,3] = r.n2 + .01807*harmonic + harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/gmr.n2)
  Zcond[2,3] = Zcond[3,2] = Zcond[1,3] = Zcond[3,1] = .01807*harmonic + harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/d.ab)
} else {
Zcond = array(0i,c(2,2)) # ohms per 1000 feet
}
Zcond[1,1] = r + .01807*harmonic + harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/gmr)
Zcond[2,2] = r.n + .01807*harmonic + harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/gmr.n)
Zcond[1,2] = Zcond[2,1] = .01807*harmonic + harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/d.ab)
Zcond = Zcond * sectionLength # actual ohms
numberOfConductors = NROW(Zcond)

# preallocate the Ybus
n = (numberOfSections + 1)*numberOfConductors
Y = array(0i,c(n,n))

# make the Ybus from the impedances
for (i in seq(from=1,to=n-numberOfConductors,by=numberOfConductors))
Y = Y.addline(Y,Zcond,from=i,to=i+numberOfConductors)

# add the shunt grounds
for (i in seq(from=2*numberOfConductors,to=n,by=numberOfConductors))
Y = Y.addshunt(Y,Zgrnd,busnum=i) # the last one is the grounded neutral

if (neutralConductor2 != "NONE") {
  for (i in seq(from=2*numberOfConductors,to=n,by=numberOfConductors))
    Y = Y.addline(Y,.00001,from=i-1,to=i)
}
# add the sub ground
Y = Y.addshunt(Y,Zsub,busnum=numberOfConductors*subBus) # the last one is the grounded neutral

# ground the phase wire at the sub
Y = Y.addline(Y,.00001,from=1+numberOfConductors*(subBus-1),to=numberOfConductors*subBus)

# make I
Isrc = array(0i,n)
# inject evenly distributed loads
if (loadType == "evenlyDistributed") {
  from = seq(from=1+numberOfConductors*subBus,to=n,by=numberOfConductors) # phase A
  to = seq(from=numberOfConductors*(1+subBus),to=n,by=numberOfConductors) # neutral
  Isrc[from] = UnbalanceI/(numberOfSections-subBus) # amperes
  Isrc[to] = -Isrc[from]
} else if (loadType == "lumpedAtEnd") {
# inject isolated load at the end
  to = n - numberOfConductors + 1
  from = n
  Isrc[from] = Isrc[from] + UnbalanceI # amperes
  Isrc[to] = Isrc[to] - Isrc[from]
} else if (loadType == "bothCircuits") {
  to = n - numberOfConductors + 1
  from = n
  Isrc[from] = Isrc[from] + UnbalanceI/2 # amperes
  Isrc[to] = Isrc[to] - Isrc[from]
  to = numberOfConductors
  from = 1
  Isrc[from] = Isrc[from] - UnbalanceI/2 # amperes
  Isrc[to] = Isrc[to] - Isrc[from]
} else {
  Isrc[numberOfConductors*subBus] = Isrc[numberOfConductors*subBus] + UnbalanceI # amperes
}
# Find the voltages:
V = solve(Y,Isrc)

Vn = V[seq(from=numberOfConductors,to=n,by=numberOfConductors)]
x = (0:numberOfSections - subBus + 1) * sectionLength / 5.28
x2 = (1:numberOfSections - subBus) * sectionLength / 5.28
idx.p = seq(from=1,to=n,by=numberOfConductors)
idx.n = seq(from=2,to=n,by=numberOfConductors)
Vx = rbind(V[idx.p] - V[idx.p + numberOfConductors],V[idx.n] - V[idx.n + numberOfConductors])
if (neutralConductor2 != "NONE") {
idx.n2 = seq(from=3,to=n,by=numberOfConductors)
Vx = rbind(Vx,V[idx.n2] - V[idx.n2 + numberOfConductors])
}
Ycond = solve(Zcond)
I = Ycond %*% Vx
Ip = I[1,]
In = I[numberOfConductors,]
Ig = Ig = colSums(I[1:numberOfConductors,])
cat("Substation Vn =",abs(Vn[subBus]),"V\n")
newgraph(width=6,height=3)
par(mfcol=c(1,2),mar=c(3, 3, 2, 1))
plot(x,abs(Vn),type="l",ylim=c(0,2*abs(Vn[subBus])),xlab="",ylab="")
mtext("Neutral-to-earth voltage, V",side=3,line=.7)
mtext("Distance from the sub, miles",side=1,line=-1,outer=TRUE)
plot(x,abs(Ig),type="l",col="blue",ylim=c(0,UnbalanceI),xlab="",ylab="")
mtext("Currents, A",side=3,line=.7)
lines(x,abs(In),col="red")
lines(x,abs(I[2,]),col="red") # 2nd neutral if there
lines(x,abs(Ip))
legend(.5,5,c("phase","neutrals","ground"),col=c("black","red","blue"),lty=1,bty="n",cex=.7)
HTMLon()
showgraph()