library(ggplot2)


Omega <- 7.292115e-5  
a <- 6371000         

# box
df <- data.frame(
  lat_deg = seq(-89.9, 89.9, by = 0.1),  
  u_ms = numeric(1799)
)

# calc u
for (i in 1:nrow(df)) {
  phi <- df$lat_deg[i] * pi/180  # rad
  df$u_ms[i] <- (Omega * a * sin(phi)^2) / cos(phi)
}

# draw ggplot
p <- ggplot(df, aes(x = lat_deg, y = u_ms)) +
  geom_line(color = "steelblue", linewidth = 1) +
  labs(
    title = expression(paste( u == frac(Omega~a~sin^2*phi, cos*phi))),
    x = "latitude (°)",
    y = expression(paste("Zonal Wind Speed, unit:m·s"^{-1}))
  ) +
  scale_x_continuous(
    breaks = seq(-90, 90, 15),
    limits = c(-90, 90),
    expand = c(0, 0)
  ) +
  scale_y_continuous(
    breaks = seq(0, 2000, 200),
    limits = c(0, 2000),
    expand = c(0, 0)
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12),
    axis.title.y = element_text(margin = margin(r = 10)),
    axis.text = element_text(size = 10, color = "black"),
    axis.text.x.minor = element_blank(),
    panel.grid.major = element_line(color = "#E6E6E6", size = 0.3),
    panel.grid.minor = element_line(color = "#F0F0F0", size = 0.2),
    panel.background = element_rect(fill = "white", color = "black", size = 1.2),  
    plot.background = element_rect(fill = "white", color = NA),
    axis.line = element_line(color = "black", size = 0.5),
    axis.ticks = element_line(color = "black"),
    axis.ticks.length = unit(0.2, "cm")
  ) +

  geom_vline(xintercept = 0, linetype = "dashed", color = "darkred", alpha = 0.7) +
  annotate("text", x = 2, y = 1900, label = "Equator", color = "darkred", hjust = 0)

print(p)

ggsave("zonalu.pdf", plot = p, device = cairo_pdf, 
       width = 8, height = 6, units = "in", dpi = 300)