"""
km_project_coords.py
For paths where B-side has only integer-degree coordinates, compute the
B-site location by projecting from the A-site using the bearing and path
length filed in PA.dat. Patches the network definition CSV.

Run after build + cleanup:
    python km_project_coords.py
"""
import os, math
import pandas as pd

ULS_DIR  = r"D:\FCC_ULS"
CSV      = r"D:\KinderMorgan\Kinder_Morgan_Network_Definition_v1.csv"
THRESHOLD = 55.0   # paths over this are suspect

def haversine_miles(lat1,lon1,lat2,lon2):
    R=3958.8
    lat1,lon1,lat2,lon2=map(math.radians,[lat1,lon1,lat2,lon2])
    dlat=lat2-lat1; dlon=lon2-lon1
    a=math.sin(dlat/2)**2+math.cos(lat1)*math.cos(lat2)*math.sin(dlon/2)**2
    return R*2*math.asin(math.sqrt(a))

def calc_bearing(lat1,lon1,lat2,lon2):
    lat1,lon1,lat2,lon2=map(math.radians,[lat1,lon1,lat2,lon2])
    dlon=lon2-lon1
    x=math.sin(dlon)*math.cos(lat2)
    y=math.cos(lat1)*math.sin(lat2)-math.sin(lat1)*math.cos(lat2)*math.cos(dlon)
    return (math.degrees(math.atan2(x,y))+360)%360

def project_point(lat, lon, bearing_deg, distance_miles):
    """Given a start point, bearing (degrees true), and distance (miles),
    compute the destination lat/lon."""
    R = 3958.8
    d = distance_miles / R
    lat1 = math.radians(lat)
    lon1 = math.radians(lon)
    brg  = math.radians(bearing_deg)
    lat2 = math.asin(math.sin(lat1)*math.cos(d) +
                     math.cos(lat1)*math.sin(d)*math.cos(brg))
    lon2 = lon1 + math.atan2(math.sin(brg)*math.sin(d)*math.cos(lat1),
                              math.cos(d)-math.sin(lat1)*math.sin(lat2))
    return math.degrees(lat2), math.degrees(lon2)

# ── Load PA to get filed bearing and path length ──────────────────────────────
print("Loading PA.dat...")
pa = pd.read_csv(os.path.join(ULS_DIR,"PA.dat"), sep="|", header=None,
                 dtype=str, on_bad_lines="skip", encoding="latin-1")
# col[4]=tx_cs col[6]=tx_loc col[7]=tx_ant col[16]=rx_cs
# col[17]=bearing col[18]=path_length (in miles per FCC spec)
pa_cs  = pa.iloc[:,4].str.strip()
pa_loc = pa.iloc[:,6].str.strip()
pa_ant = pa.iloc[:,7].str.strip()
pa_rxcs= pa.iloc[:,16].str.strip() if pa.shape[1]>16 else ""
pa_brg = pd.to_numeric(pa.iloc[:,17] if pa.shape[1]>17 else None, errors="coerce")
pa_len = pd.to_numeric(pa.iloc[:,18] if pa.shape[1]>18 else None, errors="coerce")

pa_df = pd.DataFrame({
    "tx_cs": pa_cs, "tx_loc": pa_loc, "tx_ant": pa_ant,
    "rx_cs": pa_rxcs, "bearing": pa_brg, "path_length": pa_len
})

# Build lookup: (tx_cs, rx_cs) -> (bearing, path_length)
pa_df["key"] = pa_df["tx_cs"] + "|" + pa_df["rx_cs"]
pa_lookup = pa_df[pa_df["bearing"].notna() & pa_df["path_length"].notna()]
pa_lookup = pa_lookup.set_index("key")[["bearing","path_length"]]
print(f"  PA records with bearing+length: {len(pa_lookup)}")

# Sample to verify units
sample = pa_lookup.head(10)
print(f"  Sample bearings: {sample['bearing'].tolist()}")
print(f"  Sample path lengths: {sample['path_length'].tolist()}")

# ── Load and patch CSV ────────────────────────────────────────────────────────
df = pd.read_csv(CSV, dtype=str)
df["Path_Length_Miles"] = pd.to_numeric(df["Path_Length_Miles"], errors="coerce")
df["A_Latitude"]  = pd.to_numeric(df["A_Latitude"],  errors="coerce")
df["A_Longitude"] = pd.to_numeric(df["A_Longitude"], errors="coerce")
df["B_Latitude"]  = pd.to_numeric(df["B_Latitude"],  errors="coerce")
df["B_Longitude"] = pd.to_numeric(df["B_Longitude"], errors="coerce")
df["A_Bearing"]   = pd.to_numeric(df["A_Bearing"],   errors="coerce")

patched = 0
for idx, row in df.iterrows():
    if row["Path_Length_Miles"] <= THRESHOLD:
        continue

    a_cs = row["A_Call_Sign"]; b_cs = row["B_Call_Sign"]
    a_lat = row["A_Latitude"]; a_lon = row["A_Longitude"]

    if pd.isna(a_lat) or pd.isna(a_lon):
        continue

    # Look up filed bearing and path length from PA
    key = f"{a_cs}|{b_cs}"
    if key not in pa_lookup.index:
        # Try reverse direction
        key2 = f"{b_cs}|{a_cs}"
        if key2 in pa_lookup.index:
            # reverse: bearing from B to A, so A-to-B bearing is opposite
            filed_brg = (pa_lookup.loc[key2,"bearing"] + 180) % 360
            filed_len = pa_lookup.loc[key2,"path_length"]
        else:
            print(f"  {row['path_num']} {a_cs}->{b_cs}: no PA bearing found")
            continue
    else:
        filed_brg = pa_lookup.loc[key,"bearing"]
        filed_len = pa_lookup.loc[key,"path_length"]

    if pd.isna(filed_brg) or pd.isna(filed_len) or filed_len <= 0:
        print(f"  {row['path_num']} {a_cs}->{b_cs}: invalid bearing={filed_brg} len={filed_len}")
        continue

    # PA path length is in km per FCC spec -- convert to miles
    filed_len_miles = filed_len * 0.621371

    # Project B-site from A-site
    b_lat_new, b_lon_new = project_point(a_lat, a_lon, filed_brg, filed_len_miles)
    old_mi = row["Path_Length_Miles"]
    new_mi = round(haversine_miles(a_lat, a_lon, b_lat_new, b_lon_new), 3)

    df.at[idx, "B_Latitude"]       = round(b_lat_new, 6)
    df.at[idx, "B_Longitude"]      = round(b_lon_new, 6)
    df.at[idx, "Path_Length_Miles"] = new_mi
    df.at[idx, "A_Bearing"]        = round(filed_brg, 2)
    df.at[idx, "B_Bearing"]        = round((filed_brg + 180) % 360, 2)
    print(f"  {row['path_num']} {a_cs}->{b_cs}: {old_mi:.1f} -> {new_mi:.1f} mi  brg={filed_brg:.1f}°")
    patched += 1

print(f"\nPatched {patched} paths")
df.to_csv(CSV, index=False)
print(f"Saved -> {CSV}")
print(f"\nPath length distribution after patch:")
print(df["Path_Length_Miles"].describe())
print(f"\nPaths still over {THRESHOLD} miles:")
still = df[df["Path_Length_Miles"] > THRESHOLD]
if still.empty:
    print("  None.")
else:
    print(still[["path_num","A_Call_Sign","B_Call_Sign","Path_Length_Miles"]].to_string())