#!/usr/bin/env python3
"""Extract force coefficients from Phase 3 yaw sweep.
Uses raw forces from forces1 and decomposes into drag/side force
relative to the wind direction.
"""
import os, glob, math

basedir = os.path.dirname(os.path.abspath(__file__))
angles = [0, 5, 10, 15, 20]

# Reference values
rho = 1.225      # kg/m^3
U = 10.0         # m/s
Aref = 0.006648  # m^2 (frontal area = pi * 0.046^2)
q = 0.5 * rho * U**2  # dynamic pressure

results = []
print(f"{'Yaw':>4s}  {'Fx(N)':>10s}  {'Fy(N)':>10s}  {'Fz(N)':>10s}  "
      f"{'Drag(N)':>10s}  {'Side(N)':>10s}  {'Cd':>8s}  {'Cs':>8s}  "
      f"{'Mz(Nm)':>10s}  {'CmYaw':>8s}")

for angle in angles:
    # Read forces from forces1
    forces_dir = os.path.join(basedir, f"postProcessing_yaw{angle}", "forces1", "0")
    force_file = os.path.join(forces_dir, "force.dat")
    if not os.path.isfile(force_file):
        # Try alternate names
        candidates = glob.glob(os.path.join(forces_dir, "*.dat"))
        force_file = candidates[0] if candidates else None

    moment_file = os.path.join(forces_dir, "moment.dat")
    if not os.path.isfile(moment_file):
        candidates = glob.glob(os.path.join(forces_dir, "moment*.dat"))
        moment_file = candidates[0] if candidates else None

    # Also try the coefficient file for Cd directly
    coeff_file = os.path.join(basedir, f"postProcessing_yaw{angle}",
                              "forceCoeffs1", "0", "coefficient.dat")

    # Parse forces - format: Time ((px py pz) (vx vy vz) (px py pz))
    # Actually in OF2406, force.dat has: Time total_x total_y total_z pressure_x ...
    if force_file and os.path.isfile(force_file):
        last_line = None
        with open(force_file) as f:
            for line in f:
                if not line.startswith('#'):
                    last_line = line.strip()

        if last_line:
            # Parse - the format might have parentheses
            # Remove parentheses and split
            clean = last_line.replace('(', '').replace(')', '').split()
            # Determine format by checking header
            with open(force_file) as f:
                header = ""
                for line in f:
                    if line.startswith('#'):
                        header = line
                    else:
                        break

            # Try to get total forces
            # Format: Time Fx Fy Fz ...  or Time ((Fx Fy Fz) (Fx Fy Fz) ...)
            vals = clean
            Fx_total = float(vals[1])
            Fy_total = float(vals[2])
            Fz_total = float(vals[3])
        else:
            Fx_total = Fy_total = Fz_total = float('nan')
    else:
        Fx_total = Fy_total = Fz_total = float('nan')

    # Parse moments
    Mz_total = float('nan')
    if moment_file and os.path.isfile(moment_file):
        last_line = None
        with open(moment_file) as f:
            for line in f:
                if not line.startswith('#'):
                    last_line = line.strip()
        if last_line:
            clean = last_line.replace('(', '').replace(')', '').split()
            Mz_total = float(clean[3])  # Mz

    # Decompose into drag (along wind) and side force (perpendicular)
    yaw_rad = math.radians(angle)
    cos_a = math.cos(yaw_rad)
    sin_a = math.sin(yaw_rad)

    # Wind direction = (cos_a, sin_a, 0)
    # Drag = F . wind_dir (force component along wind)
    Fdrag = Fx_total * cos_a + Fy_total * sin_a
    # Side force = F . perp_dir, where perp = (-sin_a, cos_a, 0)
    Fside = -Fx_total * sin_a + Fy_total * cos_a

    Cd = Fdrag / (q * Aref)
    Cs = Fside / (q * Aref)
    CmYaw = Mz_total / (q * Aref * 0.092) if not math.isnan(Mz_total) else float('nan')

    results.append({
        'yaw': angle, 'Fx': Fx_total, 'Fy': Fy_total, 'Fz': Fz_total,
        'Fdrag': Fdrag, 'Fside': Fside,
        'Cd': Cd, 'Cs': Cs, 'Mz': Mz_total, 'CmYaw': CmYaw
    })

    print(f"{angle:4d}  {Fx_total:10.4f}  {Fy_total:10.4f}  {Fz_total:10.4f}  "
          f"{Fdrag:10.4f}  {Fside:10.4f}  {Cd:8.4f}  {Cs:8.4f}  "
          f"{Mz_total:10.6f}  {CmYaw:8.4f}")

# Write CSV
csv_path = os.path.join(basedir, "results.csv")
with open(csv_path, 'w') as f:
    f.write("yaw_deg,Fx_N,Fy_N,Fz_N,Fdrag_N,Fside_N,Cd,Cs,Mz_Nm,CmYaw\n")
    for r in results:
        f.write(f"{r['yaw']},{r['Fx']:.6e},{r['Fy']:.6e},{r['Fz']:.6e},"
                f"{r['Fdrag']:.6e},{r['Fside']:.6e},"
                f"{r['Cd']:.6f},{r['Cs']:.6f},{r['Mz']:.6e},{r['CmYaw']:.6f}\n")
print(f"\nResults written to {csv_path}")