#----------------------------------------------------    
#  Forza di una carica allineata due cariche AVENTI STESSO SEGNO
#  in funzione della posizione
#    
#  GdA 9/10/2021
#----------------------------------------------------

q  <- -1   
Q1 <- 1
Q2 <- 4
# se segni opposti => entrambi positivi!
if(Q1*Q2 < 0) { Q1 <- abs(Q1); Q2 <- abs(Q2) } 

alpha <- sqrt(Q1/Q2) 

cat(sprintf("Q1 = %.3f,  Q2 = %.3f \n", Q1, Q2))
cat(sprintf("alpha = %.3f \n", alpha))

D <- 1    #  posizione di Q2, mentre Q1 è posizionata in x=0

k <- 1    # per semplicità usiamo unità arbitrarie (idem per le distanze)

x.eq <- D / (1 + 1/alpha)
x.sp <- D / (1 - 1/alpha)

# stampe dei risultati --------------------------------------------------
cat(sprintf("Punto equilibrio: x = %.2f\n", x.eq))
cat(sprintf("Soluzione spuria: x = %.2f\n\n", x.sp))
cat(sprintf("F1(equil.): = %.3f;  F2(equil.): = %.3f \n",
            k*q*Q1/abs(x.eq)^3*x.eq,  - k*q*Q2/abs(D-x.eq)^3*(D-x.eq)))
cat(sprintf("F1(spuria): = %.3f;  F2(spuria): = %.3f \n",
            k*q*Q1/abs(x.sp)^3*x.sp,  - k*q*Q2/abs(D-x.sp)^3*(D-x.sp)))    

# grafica ----------------------------------------------------------------
F.max <- 20 # parametro per la scala verticale
curve(k*q*Q1/abs(x)^3*x -k*q*Q2/abs(D-x)^3*(D-x), 0, D, lwd=2, 
      xlim=c(0, D), col='blue', ylab='F', ylim=c(-F.max,F.max))
curve(k*q*Q1/abs(x)^3*x, 0, D, col='orange', lwd=2, add=TRUE)
curve(-k*q*Q2/abs(D-x)^3*(D-x), -0.2, D, col='cyan', lwd=2, add=TRUE)
points(0, 0, pch=19, cex=2, col=ifelse(Q1>0, "red", "black"))
points(D, 0, pch=19, cex=2, col=ifelse(Q2>0, "red", "black"))
points(x.eq, 0, pch=19, cex=1.3, col=ifelse(q>0, "red", "black"))
abline(h=0)
abline(v=x.eq, col='green', lty=2)