library(shiny)
library(ggplot2)
library(grid)

# ---- Pearson Type III standardized quantile (negative-skew bug fixed) ----
pearson3_std_quantile <- function(p, g) {
  if (abs(g) < 1e-10) return(qnorm(p))
  alpha <- 4 / g^2
  p_adj <- if (g > 0) p else 1 - p          # flip lookup direction for g < 0
  y     <- qgamma(p_adj, shape = alpha, scale = 1)
  sign(g) * (y - alpha) / sqrt(alpha)
}

pearson3_quantile <- function(p, meanlog, sdlog, skew) {
  meanlog + sdlog * pearson3_std_quantile(p, skew)
}

aep_to_z <- function(aep) qnorm(1 - aep)

# ---- Reference return periods (computed once at startup) ----
rp_ref <- c(1.11111, 2, 5, 10, 100, 500, 1000)
rp_z   <- aep_to_z(1 / rp_ref)
rp_lbl <- c("1.1-yr", "2-yr", "5-yr", "10-yr", "100-yr", "500-yr", "1,000-yr")

# ---- UI ----
ui <- fluidPage(
  tags$head(tags$style(HTML("
    body { background: #ffffff; }
    .panel {
      background: #fafafa; border: 1px solid #e6e6e6; border-radius: 10px;
      padding: 16px 18px; margin-bottom: 14px;
    }
    .panel-base    { border-left: 5px solid #999; }
    .panel-current { border-left: 5px solid #1f4e79; }
    .curve-label {
      font-weight: 700; font-size: 14px;
      margin-bottom: 10px; padding-bottom: 8px;
      border-bottom: 1px solid #e0e0e0;
    }
    .note  { color: #555; font-size: 12px; line-height: 1.4; margin-bottom: 10px; }
    .title { font-weight: 700; font-size: 20px; margin: 6px 0 16px; }
  "))),
  
  div(class = "title", "Log-Pearson Type III Frequency Curve"),
  
  fluidRow(
    
    # ---- Left column: controls ----
    column(4,
           div(class = "note",
               "Typical ranges \u2014 ",
               "Mean: 2\u20135 (log\u2081\u2080 cfs) \u00b7 ",
               "Std dev: 0.10\u20130.50 \u00b7 ",
               "Skew: \u22120.5 to +1.0."
           ),
           
           # Baseline panel
           div(class = "panel panel-base",
               div(class = "curve-label", style = "color: #888;",
                   "\u2015\u2015  Baseline curve  (grey)"),
               numericInput("base_meanlog",
                            HTML("Mean of log<sub>10</sub>(Q)"),
                            value = 3.0, step = 0.05),
               numericInput("base_sdlog",
                            HTML("Std dev of log<sub>10</sub>(Q)"),
                            value = 0.30, min = 0.01, step = 0.01),
               numericInput("base_skew",
                            HTML("Skew of log<sub>10</sub>(Q)"),
                            value = 0.20, step = 0.05)
           ),
           
           # Current panel
           div(class = "panel panel-current",
               div(class = "curve-label", style = "color: #1f4e79;",
                   "\u2015\u2015  Current curve  (blue)"),
               numericInput("meanlog",
                            HTML("Mean of log<sub>10</sub>(Q)"),
                            value = 3.0, step = 0.05),
               numericInput("sdlog",
                            HTML("Std dev of log<sub>10</sub>(Q)"),
                            value = 0.30, min = 0.01, step = 0.01),
               numericInput("skew",
                            HTML("Skew of log<sub>10</sub>(Q)"),
                            value = 0.20, step = 0.05)
           )
    ),
    
    # ---- Right column: plot ----
    column(8,
           div(class = "panel",
               uiOutput("flow_summary"),
               plotOutput("p", height = "520px")
           )
    )
  )
)

# ---- Server ----
server <- function(input, output, session) {
  
  make_curve <- function(meanlog, sdlog, skew) {
    aep  <- seq(0.001, 0.999, length.out = 900)
    xlog <- pearson3_quantile(1 - aep, meanlog, sdlog, skew)
    data.frame(aep = aep, z = aep_to_z(aep), q = 10^xlog)
  }
  
  baseline_df <- reactive({
    make_curve(input$base_meanlog, input$base_sdlog, input$base_skew)
  })
  current_df <- reactive({
    make_curve(input$meanlog, input$sdlog, input$skew)
  })
  
  # ---- Replace the aep_ticks block and plot internals only ----
  
  output$flow_summary <- renderUI({
    fmt <- function(q) formatC(round(q), format = "d", big.mark = ",")
    
    label_aeps <- c(1/100, 1/1000)
    bq <- 10^pearson3_quantile(1 - label_aeps,
                               input$base_meanlog, input$base_sdlog, input$base_skew)
    cq <- 10^pearson3_quantile(1 - label_aeps,
                               input$meanlog,      input$sdlog,      input$skew)
    
    tags$div(
      style = "display:flex; gap:32px; margin-bottom:10px;
               font-size:13px; font-family:monospace;",
      tags$div(
        style = "color:#888; border-left:3px solid #bbb; padding-left:8px;",
        tags$div(tags$b("Baseline")),
        tags$div(paste0("100-yr:     ", fmt(bq[1]), " cfs")),
        tags$div(paste0("1,000-yr:  ", fmt(bq[2]), " cfs"))
      ),
      tags$div(
        style = "color:#1f4e79; border-left:3px solid #1f4e79; padding-left:8px;",
        tags$div(tags$b("Current")),
        tags$div(paste0("100-yr:     ", fmt(cq[1]), " cfs")),
        tags$div(paste0("1,000-yr:  ", fmt(cq[2]), " cfs"))
      )
    )
  })
  
  output$p <- renderPlot({
    base <- baseline_df()
    cur  <- current_df()
    
    # ---- Axis ticks ----
    rp_ticks   <- c(1.111, 2, 5, 10, 100, 500, 1000)
    aep_ticks  <- 1 / rp_ticks
    z_breaks   <- aep_to_z(aep_ticks)
    aep_labels <- c("0.9", "0.5", "0.2", "0.1", "0.01", "0.002", "0.001")
    rp_labels  <- c("1.1", "2", "5", "10", "100", "500", "1,000")
    
    # ---- Flow estimates at 100-yr and 1,000-yr ----
    label_aeps <- c(1/100, 1/1000)
    fmt <- function(q) formatC(round(q), format = "d", big.mark = ",")
    bq  <- 10^pearson3_quantile(1 - label_aeps,
                                input$base_meanlog, input$base_sdlog, input$base_skew)
    cq  <- 10^pearson3_quantile(1 - label_aeps,
                                input$meanlog,      input$sdlog,      input$skew)
    
    
    # ---- Points on curve at 100-yr and 1,000-yr ----
    label_aeps <- c(1/100, 1/1000)
    label_z    <- aep_to_z(label_aeps)
    
    base_pts <- data.frame(z = label_z, q = bq)
    cur_pts  <- data.frame(z = label_z, q = cq)
    
    ggplot() +
      geom_line(data = base, aes(x = z, y = q),
                linewidth = 1.0, color = "grey65") +
      geom_line(data = cur,  aes(x = z, y = q),
                linewidth = 1.2, color = "#1f4e79") +
      geom_point(data = base_pts, aes(x = z, y = q),
                 color = "grey55", size = 3, shape = 19) +
      geom_point(data = cur_pts,  aes(x = z, y = q),
                 color = "#1f4e79", size = 3, shape = 19) +
      scale_y_log10() +
      scale_x_continuous(
        name   = "Annual Exceedance Probability (AEP)",
        breaks = z_breaks,
        labels = aep_labels,
        sec.axis = sec_axis(
          trans  = ~ .,
          name   = "Return Period (years)",
          breaks = z_breaks,
          labels = rp_labels
        )
      ) +
      labs(y = "Peak Flow (cfs) [log scale]") +
      theme_minimal(base_size = 15) +
      theme(
        panel.grid.minor    = element_blank(),
        axis.title.x.bottom = element_text(margin = margin(t = 8)),
        axis.title.x.top    = element_text(margin = margin(b = 8))
      )
  })
}

shinyApp(ui, server)