#!/usr/bin/env bash
set -euo pipefail

usage() {
    cat <<'EOF'
Usage:
  ./run_ff_slip_dist_latlon.sh one
  ./run_ff_slip_dist_latlon.sh two
  ./run_ff_slip_dist_latlon.sh both

Optional environment overrides:
  ONE_DATA_DIR=/path/to/one_segment_inversion_directory
  TWO_DATA_DIR=/path/to/two_segment_inversion_directory
  PYTHON_CMD=python3
  PLOT_SCRIPT=/path/to/ff_slip_dist_latlon.py
  ONE_COMBINED_FILE=/path/to/one_segment.plt
  TWO_COMBINED_FILE=/path/to/two_segment.plt
  ONE_FRW_FILE=/path/to/one_segment.frw
  TWO_FRW_FILE=/path/to/two_segment.frw
  NUM_TIME_WINDOWS=3
  REQUIRE_SLAB_PROFILE=0  # 1 makes missing Slab2 coverage fatal
  MAKE_FRW_PATCHES=1      # default: generate configured patches and show boundaries
EOF
}

[[ $# -eq 1 ]] || { usage; exit 2; }
CASE_NAME=$1

#PYTHON_CMD=${PYTHON_CMD:-python}
PYTHON_CMD=${PYTHON_CMD:-/data/seis01/taira/miniconda3/envs/netops/bin/python}

PLOT_SCRIPT=${PLOT_SCRIPT:-ff_slip_dist_latlon.py}

# Separate directories containing the one- and two-segment reg_inv .plt files.
# Override either directory independently for another inversion result.
ONE_DATA_DIR=${ONE_DATA_DIR:-/data/seis02/taira/BSL_MT_Training/mendo2024/ff_loop_test_rakeval/work_ff_inc_gps_rakeval_multime3.loop_strike_loop/inv_para_277.00_84.00_-135.00+-225.00_40_20_80.00_40.00_40.374_-125.022_9.0_5.00_4.00_0_0_2_567/source_para_2.80_1.00_0.0002.rakeval_multime3}
TWO_DATA_DIR=${TWO_DATA_DIR:-/data/seis02/taira/BSL_MT_Training/mendo2024/ff_loop_test_rakeval/work_ff_inc_gps_seg2_rdf3_side3_rakeval_multime3_loop/inv_para_275.00_84.00_-135.00+-225.00_20_20_40.00_40.00_40.374_-125.022_9.0_-5.00_4.00_0_0_2_567/source_para_2.20_1.00_0.0002.seg2_rdf3_side3_rakeval_multime3/280.00_84.00_-135.00+-225.00_20_20_40.00_40.00_40.330560_-124.695681_1.043825_10.00_4.00_0/source_para2_2.80_delays_0.0_0.0_0.5}

# Each full input path can still be overridden independently.
ONE_COMBINED_FILE=${ONE_COMBINED_FILE:-${ONE_DATA_DIR}/reg_inv.out.inc_gps_rakeval_multime30.plt}
TWO_COMBINED_FILE=${TWO_COMBINED_FILE:-${TWO_DATA_DIR}/reg_inv.out.inc_gps_seg2_rdf3_side3_rakeval_multime30.plt}
ONE_FRW_FILE=${ONE_FRW_FILE:-${ONE_DATA_DIR}/reg_inv.out.inc_gps_rakeval_multime3.frw}
TWO_FRW_FILE=${TWO_FRW_FILE:-${TWO_DATA_DIR}/reg_inv.out.inc_gps_seg2_rdf3_side3_rakeval_multime3.frw}
NUM_TIME_WINDOWS=${NUM_TIME_WINDOWS:-3}
REQUIRE_SLAB_PROFILE=${REQUIRE_SLAB_PROFILE:-0}
MAKE_FRW_PATCHES=${MAKE_FRW_PATCHES:-1}
case "$REQUIRE_SLAB_PROFILE" in
    0|1) ;;
    *) echo "ERROR: REQUIRE_SLAB_PROFILE must be 0 or 1" >&2; exit 2 ;;
esac
case "$MAKE_FRW_PATCHES" in
    0|1) ;;
    *) echo "ERROR: MAKE_FRW_PATCHES must be 0 or 1" >&2; exit 2 ;;
esac

# Edit these arrays to define one, two, or three independent forward-model patches.
# Format: --frw_patchN SEGMENT XMIN_KM XMAX_KM ZMIN_KM ZMAX_KM
# The one-segment case can only use segment ID 1.
ONE_FRW_PATCH_ARGS=(
    --frw_patch1 1 56 74 0 20
    --frw_patch1_output one_fault_patch1.frw
    --frw_patch2 1 16 34 22 32
    --frw_patch2_output one_fault_patch2.frw
    --frw_patch3 1 38 50 6 14
    --frw_patch3_output one_fault_patch3.frw
    --frw_patch_summary one_fault_patch_summary.csv
    --frw_patch_beachball_file one_frw_patch_focal_mechanisms.pdf
)

TWO_FRW_PATCH_ARGS=(
    --frw_patch1 1 13 35 -1 21
    --frw_patch1_output fault_patch1.frw
    --frw_patch2 2 7 35 9 25
    --frw_patch2_output fault_patch2.frw
    --frw_patch3 1 -1 11 5 15
    --frw_patch3_output fault_patch3.frw
    --frw_patch_summary fault_patch_summary.csv
    --frw_patch_beachball_file two_frw_patch_focal_mechanisms.pdf
)

SEIS_CATALOG=../hypodd_catalog/NCAeqDDRT.v202201.downloaded.20260624
REPEATING_CATALOG=MTJ.freq4.0-16.0Hz_maxdist60.0_coh9700_linkage_cluster.txt
SLAB_DIR=/data/seis02/taira/BSL_MT_Training/mendo2024/ff_plot/Slab2Distribute_Mar2018
GMT_CMD=${GMT_CMD:-gmt}

require_file() {
    local label=$1
    local path=$2
    [[ -f "$path" ]] || { echo "ERROR: $label not found: $path" >&2; exit 1; }
}

require_dir() {
    local label=$1
    local path=$2
    [[ -d "$path" ]] || { echo "ERROR: $label directory not found: $path" >&2; exit 1; }
}

preflight_common() {
    [[ -x "$PYTHON_CMD" ]] || { echo "ERROR: Python executable not found or not executable: $PYTHON_CMD" >&2; exit 1; }
    require_file "plot script" "$PLOT_SCRIPT"
    require_file "seismicity catalog" "$SEIS_CATALOG"
    require_file "repeating-earthquake catalog" "$REPEATING_CATALOG"
    require_dir "Slab2" "$SLAB_DIR"
    command -v "$GMT_CMD" >/dev/null 2>&1 || { echo "ERROR: GMT command not found: $GMT_CMD" >&2; exit 1; }

    local script_dir
    script_dir=$(cd "$(dirname "$PLOT_SCRIPT")" && pwd)
    PYTHONPATH="$script_dir${PYTHONPATH:+:$PYTHONPATH}" "$PYTHON_CMD" -c "import ff_catalog_common, ff_map_common" \
        || { echo "ERROR: could not import ff_catalog_common and ff_map_common from $script_dir or PYTHONPATH" >&2; exit 1; }
}

COMMON_ARGS=(
    --make_plot
    --use_moment_for_slip
    --slip_cmax 400
    --marker_size 125
    --plot_mu_depth
    --plot_mu_compare

    --plot_seismicity
    --seis_catalog "$SEIS_CATALOG"
    --seis_start_time 2024-12-05T18:44:21
    --seis_end_time 2025-12-05T18:44:21
    --seis_min_depth 0
    --seis_max_depth 50
    --seis_min_mag 1.5
    --seis_edgecolor blue
    --seis_facecolor none
    --seis_max_across_km 20
    # Keep every event on every segment whose +/-20 km corridor test passes.
    # The same earthquake may therefore be projected on both segment panels.
    --seis_across_filter
    --seis_segment_assignment_mode all_candidates
    --seis_marker_ddrt30 o
    --seis_marker_hinv D
    --seis_marker_other x
    --seis_scale_by_mag
    --seis_marker_size 30
    --seis_marker_mag_ref 3.0
    --seis_marker_mag_scale 3.0
    --seis_marker_size_min 8
    --seis_marker_size_max 220
    --seis_depth_offset_km 0.293
    # Explicitly use the same datum correction for Slab2; do not inherit it silently.
    --slab_depth_offset_km 0.293
    --ignore_across_strike_distance_depth_projection
    --seis_legend
    --seis_legend_loc "lower left"

    --plot_repeating
    --repeating_catalog "$REPEATING_CATALOG"
    --repeating_facecolor red

    --plot_slab
    --slabdir "$SLAB_DIR"
    --slab_id cas
    --slab_color black
    --slab_sample_backend gmt
    --gmt_cmd "$GMT_CMD"
    --slab_sample_dx_km 1.0

    --make_map_plot
    --map_plot_seismicity
    --map_plot_repeating

    # Show the exact +/- across-fault strips used to select cross-section events.
    --map_plot_selection_polygons
    --map_selection_segment1_color blue
    --map_selection_segment2_color orange
    --map_selection_polygon_alpha 0.12
    --map_selection_polygon_linewidth 1.2
    --map_selection_polygon_linestyle dashed
    --map_selection_draw_center

    --no_map_plot_slip
    #--map_plot_slip


    --map_plot_target_event \
    --target_event_lat 40.374 \
    --target_event_lon -125.022 \
    --target_event_model_depth_km 9.0 \
    --target_event_depth_mode model_z \
    
    --target_event_cross_section

    # Fixed bounds for consistent one/two-segment comparisons
    # --map_min_lon -125.4
    # --map_max_lon -124.5
    # --map_min_lat 40.30
    # --map_max_lat 40.45

    # --map_min_lon -125.5
    # --map_max_lon -123.5
    # --map_min_lat 40.00
    # --map_max_lat 41.00

    --map_min_lon -125.6
    --map_max_lon -124.2
    --map_min_lat 40.10
    --map_max_lat 40.50

    # Slip cells reconstructed from actual inversion lon/lat coordinates
    --map_fault_render midpoint_cells
    --map_midpoint_cell_edgecolor none
    --map_midpoint_cell_linewidth 0
    --map_midpoint_cell_alpha 1.0
    --map_midpoint_cell_zorder 4

    # Segment outlines
    --map_midpoint_draw_outlines
    --map_segment1_outline_color black
    --map_segment2_outline_color gray
    --map_segment1_outline_linewidth 1.25
    --map_segment2_outline_linewidth 1.25
    --map_segment1_outline_linestyle solid
    --map_segment2_outline_linestyle dashed

    --map_segment_outline_zorder 14

    # Top trace is unnecessary when full outlines are shown
    --no_map_draw_top_trace

    # Map scale
    --map_scale_bar_km 20
    --map_scale_bar_loc lower_right
    --map_scale_bar_x 0.08
    --map_scale_bar_y 0.06


    # Common absolute marker areas
    --seis_legend_marker_size 48
    --map_legend_marker_size 48

    # Keep the older global multipliers neutral
    --seis_legend_markerscale 1.0
    --map_legend_markerscale 1.0

    # Per-symbol visual compensation
    --legend_scale_hinv 1.0
    --legend_scale_ddrt30 1.0
    --legend_scale_repeating_hinv 1.0
    --legend_scale_repeating_ddrt30 1.0
    --legend_scale_target 2.0


    # Coastline
    --map_coastline_source cartopy
    --map_coastline_resolution 10m
    --map_coastline_color 0.35
    --map_coastline_linewidth 0.8

    # Longitude/latitude ticks and guide lines
    --map_lon_tick_interval 0.2
    --map_lat_tick_interval 0.05
    --map_tick_label_size 9
    --map_gridlines
    --map_gridline_color 0.75
    --map_gridline_linewidth 0.5
    --map_gridline_alpha 0.6
    --map_gridline_linestyle dotted



    # Cross-section legend
    --cross_section_legend_mode segment1_combined

    --legend_label_hinv "hypoinv"
    --legend_label_ddrt30 "DDRT"

    --legend_label_repeating_hinv "repeating EQ (hypoinv)"
    --legend_label_repeating_ddrt30 "repeating EQ (DDRT)"

    --legend_label_slab "Slab2 interface"
    --legend_label_target "2024 M7.0 EQ"

    --map_legend_label_segment1 "segment 1"
    --map_legend_label_segment2 "segment 2"

    --map_legend_label_hinv "hypoinv"
    --map_legend_label_ddrt30 "DDRT"

    --map_legend_label_repeating_hinv "repeating EQ (hypoinv)"
    --map_legend_label_repeating_ddrt30 "repeating EQ (DDRT)"

    --map_legend_label_target "2024 M7.0 EQ"
    --map_legend_label_slab "Slab2 interface"

    --slip_contour_interval_cm 50
    --slip_contour_color 0.20
    --slip_contour_linewidth 0.75
    --slip_contour_linestyle solid
    --slip_contour_alpha 0.90
    --slip_contour_label_size 6.5

)

run_one() {
    local prefix=plot_test_one_segment
    local case_args=(
        --combined_file "$ONE_COMBINED_FILE"
        --frw_file "$ONE_FRW_FILE"
        --num_time_windows "$NUM_TIME_WINDOWS"
        --rupture_config rakeval_multime3

        --xgrid_num 40
        --zgrid_num 20
        --x_length_km 80
        --z_length_km 40

        --strike1_deg 277
        --seis_strike1_deg 277
        --dip1_deg 84

        --panel_wspace 0
        --panel_width_scale 1.0

        --output_prefix "$prefix"
        --plot_file "${prefix}.pdf"
        --seis_output "${prefix}.selected_seismicity.csv"
        --seis_candidate_output "${prefix}.seismicity_event_segment_candidates.csv"
        --repeating_output "${prefix}.selected_repeating.csv"
        --repeating_candidate_output "${prefix}.repeating_event_segment_candidates.csv"
        --slab_output "${prefix}.slab_profile.csv"

        --map_plot_file "${prefix}.map.pdf"

    )

    if [[ "$REQUIRE_SLAB_PROFILE" == 1 ]]; then
        case_args+=(--require_slab_profile)
    fi
    if [[ "$MAKE_FRW_PATCHES" == 1 ]]; then
        case_args+=("${ONE_FRW_PATCH_ARGS[@]}")
    fi

    preflight_common
    require_file "one-segment combined file" "$ONE_COMBINED_FILE"
    require_file "one-segment FRW file" "$ONE_FRW_FILE"

    echo "# Running one-segment case"
    echo "# Combined input: $ONE_COMBINED_FILE"
    echo "# FRW input:      $ONE_FRW_FILE"
    echo "# Output: $prefix"
    "$PYTHON_CMD" "$PLOT_SCRIPT" "${COMMON_ARGS[@]}" "${case_args[@]}"
}

run_two() {
    local prefix=plot_test_two_segments
    local case_args=(
        --combined_file "$TWO_COMBINED_FILE"
        --frw_file "$TWO_FRW_FILE"
        --num_time_windows "$NUM_TIME_WINDOWS"
        --rupture_config seg2_rakeval_multime3

        --xgrid_num 20
        --zgrid_num 20
        --x_length_km 40
        --z_length_km 40

        # Segment 2 inherits the segment 1 grid unless these are changed.
        --xgrid_num2 20
        --zgrid_num2 20
        --x_length_km2 40
        --z_length_km2 40

        --strike1_deg 275
        --strike2_deg 280
        --seis_strike1_deg 275
        --seis_strike2_deg 280
        --dip1_deg 84
        --dip2_deg 84

        --panel_wspace 0
        --panel_width_scale 0.75

        --output_prefix "$prefix"
        --plot_file "${prefix}.pdf"
        --seis_output "${prefix}.selected_seismicity.csv"
        --seis_candidate_output "${prefix}.seismicity_event_segment_candidates.csv"
        --repeating_output "${prefix}.selected_repeating.csv"
        --repeating_candidate_output "${prefix}.repeating_event_segment_candidates.csv"
        --slab_output "${prefix}.slab_profile.csv"

        --map_plot_file "${prefix}.map.pdf"
    )

    if [[ "$REQUIRE_SLAB_PROFILE" == 1 ]]; then
        case_args+=(--require_slab_profile)
    fi
    if [[ "$MAKE_FRW_PATCHES" == 1 ]]; then
        case_args+=("${TWO_FRW_PATCH_ARGS[@]}")
    fi

    preflight_common
    require_file "two-segment combined file" "$TWO_COMBINED_FILE"
    require_file "two-segment FRW file" "$TWO_FRW_FILE"

    echo "# Running two-segment case"
    echo "# Combined input: $TWO_COMBINED_FILE"
    echo "# FRW input:      $TWO_FRW_FILE"
    echo "# Output: $prefix"
    "$PYTHON_CMD" "$PLOT_SCRIPT" "${COMMON_ARGS[@]}" "${case_args[@]}"
}

case "$CASE_NAME" in
    one)  run_one ;;
    two)  run_two ;;
    both) run_one; run_two ;;
    *)
        echo "ERROR: unknown case: $CASE_NAME" >&2
        usage >&2
        exit 2
        ;;
esac