#-------------------------------------------------------------
# Calcolo e grafica interattiva 
# di curve di livello e gradienti
#   -- versione con stampe, anche in coordinate polari
#
# GdA 29/10/2021
#-------------------------------------------------------------

# Funzione: 
# https://web.northeastern.edu/seigen/MTH1252LevelCurves.html
myfun <- function(x,y) x*y*exp(-(x^2+y^2)/2)
# gradiente
Gr.x <- function(x,y) { y*exp(-(x^2+y^2)/2) +
                             x*y*exp(-(x^2+y^2)/2)*(-x)}
Gr.y <- function(x,y) { x*exp(-(x^2+y^2)/2) +
                            x*y*exp(-(x^2+y^2)/2)*(-y)}

# Funxione e gradiente in coordinate polari
myfun.pol <- function(r,th) r^2 * sin(th)*cos(th) * exp(-r^2/2) 
Gr.r <- function(r,th) {
    2*r * sin(th)*cos(th) * exp(-r^2/2) + r^2 * sin(th)*cos(th) * exp(-r^2/2) * (-r)    
}
Gr.th <- function(r, th) (cos(th)^2 - sin(th)^2) * r^2 * exp(-r^2/2) / r   

# Prepara i punti del 'contour plot' -- in questa versione SOLO PRIMO QUADRANTE
x <- (0:100)*pi/100
y <- (0:100)*pi/100
z <- outer(x, y, FUN='myfun')

# https://r-charts.com/correlation/contour-plot/
contour(x, y, z, col='gray', nlevel=20, lwd=1.5, labcex = 1.3, asp=TRUE)
  #      xaxs="i", yaxs="i", asp=TRUE)
text(3,3, "FINE", cex=2)
abline(h=0, lty=2)
abline(v=0, lty=2)
abline(0,1, lty=2)
points(0,0,pch=19)
cat(sprintf("Cliccare su un punto della finestra grafica (su FINE per terminare)\n"))
scala = 1
while(1) {  
  p <- locator(1)
  if((p$x[1]>2.8) & (p$y[1]> 2.8)) break
  # introduciamo px e py per semplificare
  px <- p$x[1]
  py <- p$y[1]
  cat(sprintf("\nP = (%.2f, %.2f)", px, py))
  # punto in coordinate polari --------------------
  r  <-  sqrt(px^2+py^2)
  th <-  atan2(py, px)
  th.gr <- th*180/pi
  cat(sprintf(" -->  r = %.2f, theta = %.2f (= %.1f°",
              r, th, th.gr))
  if(th > 0) {
     string =  ")\n"
  } else {
     th.gr.pos <- 360 + th.gr
     string =  sprintf(" -> %.1f°)\n", th.gr.pos) 
  }
  cat(string)
  
  # campo scalare ('V') -------------------------------------
  V <- myfun(px, py)
  cat(sprintf("Campo scalare: V = %.3f ", V))
  V.pol <- myfun.pol(r, th) 
  cat(sprintf("[= %.3f, ricalcolato dalle coord. polari]\n", V.pol))

  # gradienti
  cat(sprintf("Gradiente coord. cart.:  (%.3f, %.3f)\n", Gr.x(px,py), Gr.y(px,py)))
  # idem, dalle coordinate polari
  cat(sprintf("Gradiente coord. pol.:  (%.3f, %.3f)\n", Gr.r(r,th), Gr.th(r,th)))
  
  # grafica -----------------------------------------
  points(p$x[1], p$y[1],pch=19, cex=0.3, col='red')
  arrows(p$x[1], p$y[1],
         p$x[1] + Gr.x(p$x[1],p$y[1])*scala,
         p$y[1] + Gr.y(p$x[1],p$y[1])*scala,
         len=0.1, col='red')
  # componenti cartesiane
  arrows(p$x[1], p$y[1],
         p$x[1] + Gr.x(p$x[1],p$y[1])*scala, p$y[1], len=0.1, col='blue')
  arrows(p$x[1], p$y[1],
         p$x[1], p$y[1] + Gr.y(p$x[1],p$y[1])*scala, len=0.1, col='blue')
  # componenti polari 
  # versore lungo r
  ur = c(px, py) / r
  arrows(px, py,
         px + Gr.r(r, th)*ur[1]*scala,
         py + Gr.r(r, th)*ur[2]*scala, len=0.1, col='green')
  # versore 'lungo theta' (ortogonale a quello lungo r)
  mr = ur[2]/ur[1]
  mth = -1/mr
  uth = c(1, mth)  # non normalizzato
  uth = -1 * uth / sqrt(1+mth^2)   # il segno meno è stato messo, senza troppo conti,
                                   # facendo delle prove (c'erano due possibilità...) 
  arrows(px, py,
         px + Gr.th(r, th)*uth[1]*scala,
         py + Gr.th(r, th)*uth[2]*scala, len=0.1, col='green')
}