# =========================================================
# 🔹 CÓDIGO FINAL (OPTIMIZADO PARA RPUBS / RMARKDOWN)
# =========================================================

# 1. CARGAR DATOS
# Asegúrate de que el archivo esté en la misma carpeta que tu .Rmd
datos <- read.csv("Petroleo_Ontaro.csv", header = TRUE, sep = ";", stringsAsFactors = FALSE)

# 2. LIMPIEZA (Obligatoria)
clean_lat <- function(val) {
  val <- as.character(val)
  val <- gsub("\\.", "", val) 
  num <- as.numeric(val)
  if (is.na(num)) return(NA)
  if (num == 0) return(0)
  while (abs(num) > 90) { num <- num / 10 }
  return(num)
}

datos$SURFACE_LATITUDE_83 <- sapply(datos$SURFACE_LATITUDE_83, clean_lat)
datos$TOTAL_DEPTH <- as.numeric(gsub(",", ".", as.character(datos$TOTAL_DEPTH)))

# Filtrar datos válidos
datos_validos <- subset(datos, TOTAL_DEPTH > 0 & SURFACE_LATITUDE_83 > 40 & SURFACE_LATITUDE_83 < 60)

x <- datos_validos$TOTAL_DEPTH
y <- datos_validos$SURFACE_LATITUDE_83

# 3. MODELO MATEMÁTICO
x1 <- log(x)
y1 <- log(y)
regresionPotencial <- lm(y1 ~ x1)
b <- coef(regresionPotencial)[2]
a <- exp(coef(regresionPotencial)[1])

# 4. GRÁFICA (MÉTODO "LINES" - NO FALLA EN RPUBS) ----------------

# Ajustamos márgenes dentro del mismo bloque (sin abrir ventanas nuevas)
par(mar = c(4.5, 4.5, 3, 1))

# A) Dibujar los puntos
plot(x, y, 
     col = rgb(0, 0.5, 0.5, 0.4), 
     pch = 16, 
     main = "Regresión Potencial: Latitud vs Profundidad",
     xlab = "Profundidad (m)", 
     ylab = "Latitud (grados)",
     ylim = c(40, 48),    
     xlim = c(0, 5000),
     cex.main = 0.9,
     las = 1)

# B) Calcular la línea manualmente (Esto es lo que arregla el error)
# Creamos 500 puntos desde el mínimo hasta el máximo de X
x_linea <- seq(from = 10, to = 5000, length.out = 500)
# Calculamos la Y para esos puntos usando tu ecuación
y_linea <- a * x_linea^b

# C) Dibujar la línea sobre el gráfico
lines(x_linea, y_linea, col = "red", lwd = 3)

# D) Leyenda
legend("topright", 
       legend = paste("y =", round(a, 2), "* x ^", round(b, 5)),
       col = "red", lwd = 3, bty = "n", cex = 0.8)

# E) Rejilla
grid()

# 5. RESULTADOS
cat("\n--- ESTIMACIÓN --- \n")
## 
## --- ESTIMACIÓN ---
val_x <- 3850
pred_y <- a * val_x^b
cat(sprintf("Si la profundidad es %d m, la Latitud estimada es: %.4f°\n", val_x, pred_y))
## Si la profundidad es 3850 m, la Latitud estimada es: 42.1397°