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
}
library(RDCOMClient)
if (!exists("DSSengine")) {
  DSSengine = COMCreate("ESIDSSEngine.DSS")
  DSSengine$Start(0)
  DSSText = DSSengine[["Text"]]
  DSSText[["Command"]] = "clear"
}

load(RpadBaseFile("AA_conductors.RData"))
HTMLon()
cat("Neutral conductor: ")
HTMLselect("neutralConductor1", a$label, default=16)

numspans = 50
subBus = 1
totalLength = totalLengthMi * 5.28 # ohms/kfeet
sectionLength = totalLength/numspans # kfeet
Zgrnd = averageGround/(groundsPerKFeet * totalLength) *
        numspans # 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 feet
gmr.n = aa.n1$gmr

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

# Enter the lines in DSS:

dsstrimatrix = function(Z) 
  paste("(", Z[1,1], "|", Z[1,2], Z[2,2], ")")

"%paste%" = function(a,b) paste(a, b, "\n", sep="")
dss = "clear\n"

dss = dss %paste% "New Circuit.EX.1.2 basekv=12.47 Isc3=100000 Isc1=60000 Phases=1"
dss = dss %paste% 
  paste("New linecode.lspan nphases=2 BaseFreq=60",
        "rmatrix = ", dsstrimatrix(Re(Zcond)),
        "xmatrix = ", dsstrimatrix(Im(Zcond)),
        "cmatrix = (.1  | -0  .1)")

dss = dss %paste% 
  paste("New Line.sub", 1:numfeeders, 
        " bus1=sourcebus.1.2",
        " bus2=Bus", 1:numfeeders, "x1.1.2",
        " linecode=lspan length=1", 
        sep = "", collapse="\n")
dss = dss %paste% 
  paste("New Line.l", 1:(numspans*numfeeders), 
        " bus1=Bus", rep(1:numfeeders,each=numspans), 'x', rep(1:numspans,numfeeders), ".1.2",
        " bus2=Bus", rep(1:numfeeders,each=numspans), 'x', rep(1:numspans,numfeeders) + 1, ".1.2",
        " linecode=lspan length=1", 
        sep = "", collapse="\n")

# Enter the line grounds in DSS:

dss = dss %paste% 
  paste("New linecode.grounds nphases=1 BaseFreq=60",
        "rmatrix = ", Zgrnd,
        "xmatrix = 0. cmatrix = .1")
dss = dss %paste% 
  paste("New Line.grounds", 1:((numspans+1)*numfeeders), 
        " bus1=Bus", rep(1:numfeeders,each=numspans+1), 'x', rep(1:(numspans+1),numfeeders), ".2",
        " bus2=Bus", rep(1:numfeeders,each=numspans+1), 'x', rep(1:(numspans+1),numfeeders), ".0",
        " linecode=grounds length=1", 
        sep = "", collapse="\n")

# Enter the station ground in DSS:

dss = dss %paste% 
  paste("New linecode.subground nphases=1 BaseFreq=60",
        "rmatrix = ", Zsub,
        "xmatrix = 0. cmatrix = .1")
dss = dss %paste% 
  paste("New Line.subground",
        " bus1=sourcebus.2",
        " bus2=sourcebus.0", 
        " linecode=subground length=1", 
        sep = "")

# Enter the loads in DSS:

dss = dss %paste% 
  paste("New Load.unbalance", 1:numfeeders, 
        " bus1=Bus", 1:numfeeders, "x40.1.2",
        " Phases=1 kV=10 kVA=100 pf=1.0 Model=5",
        sep="", collapse="\n")

# Solve and export from DSS:

dss = dss %paste% "Solve"
dss = dss %paste% "Export Voltages"
dss = dss %paste% "Export Currents"

cat(file="strayvoltage.dss", dss)

DSSText[["Command"]] = "Redirect e:\\www\\Rpad_homepage_new\\strayvoltage.dss"

I = read.csv("EXP_CURRENTS.CSV", comment = "", strip.white = TRUE)
V = read.csv("EXP_VOLTAGES.CSV", comment = "", strip.white = TRUE)
idx = intersect(which(I$Terminal == 1), # AND
                grep("line.l", as.character(I$Element)))
idx = I$Terminal == 1 &
      regexpr("bus1", as.character(I$Bus)) > 0 &
      regexpr("line.grounds", as.character(I$Element)) > 0
Vn = Zgrnd * I$Magnitude[idx]
idx = I$Terminal == 1 & I$Phase == 1 &
      regexpr("bus1", as.character(I$Bus)) > 0 &
      regexpr("line.l", as.character(I$Element)) > 0 
Ip = I$Magnitude[idx] * an(I$Angle[idx])
idx = I$Terminal == 1 & I$Phase == 2 &
      regexpr("bus1", as.character(I$Bus)) > 0 &
      regexpr("line.l", as.character(I$Element)) > 0 
In = I$Magnitude[idx] * an(I$Angle[idx])
Ig = Ip + In
x = 0:(numspans-1) * sectionLength / 5.28

# Do the output:
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[-1]),
     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(Ip))
legend(.5, 5, c("phase","neutral","ground"), 
       col = c("black","red","blue"), 
       lty = 1, bty = "n", cex = .7)
HTMLon()
showgraph()