#!/usr/bin/env python3 """ Plot Cp distributions for the ONERA M6 wing at 7 spanwise sections. Compares experimental data (exp_case_2308.dat) vs LAVA unstructured results. Usage: python plot_cp_sections.py python plot_cp_sections.py --exp --lava --out """ import argparse import re from pathlib import Path import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt # --------------------------------------------------------------------------- # Section metadata: maps section index (0-based) to y/b label # --------------------------------------------------------------------------- YB_VALUES = [0.20, 0.44, 0.65, 0.80, 0.90, 0.96, 0.99] # --------------------------------------------------------------------------- # Parsers # --------------------------------------------------------------------------- def parse_lava_cutdata(filepath): """ Parse a Tecplot-style LAVA cutdata file. Returns ------- list of dict, one per zone: {"label": str, "xc": ndarray, "cp": ndarray} """ zones = [] with open(filepath) as fh: lines = fh.readlines() i = 0 while i < len(lines): line = lines[i].strip() # Detect start of a zone if line.upper().startswith("ZONE"): # Extract zone title m = re.search(r'T\s*=\s*"([^"]*)"', line) label = m.group(1) if m else f"Zone {len(zones)+1}" # Advance through zone header lines until we hit numeric data # Look for I= to know how many points to read npoints = None while i < len(lines): header_line = lines[i].strip() # Check for I= count (handles both "I=628" and "I= 34") m_i = re.search(r'\bI\s*=\s*(\d+)', header_line) if m_i: npoints = int(m_i.group(1)) # Stop searching header when we hit numeric data try: float(header_line.split()[0]) break # this line is data except (ValueError, IndexError): i += 1 if npoints is None: # Count until next ZONE or EOF as a fallback pass # Read npoints data lines xc_list, cp_list = [], [] count = 0 while i < len(lines) and (npoints is None or count < npoints): data_line = lines[i].strip() if not data_line: i += 1 continue # Stop if we hit a new header section if data_line.upper().startswith(("ZONE", "TITLE", "VARIABLES")): break parts = data_line.split() try: vals = [float(p) for p in parts] xc_list.append(vals[0]) # column 0: X/C cp_list.append(vals[6]) # column 6: PCOEF count += 1 except (ValueError, IndexError): pass i += 1 zones.append({ "label": label, "xc": np.array(xc_list), "cp": np.array(cp_list), }) else: i += 1 return zones def parse_experiment(filepath): """ Parse the ONERA M6 experimental data file (Tecplot POINT format). Variables: Section, Tap, X/L, Z/L, CP Returns ------- list of dict, one per zone: {"label": str, "xc": ndarray, "cp": ndarray} """ zones = [] with open(filepath) as fh: lines = fh.readlines() i = 0 while i < len(lines): line = lines[i].strip() if line.upper().startswith("ZONE"): m_label = re.search(r'T\s*=\s*"([^"]*)"', line) label = m_label.group(1).strip().lstrip(",").strip() if m_label else f"Zone {len(zones)+1}" m_i = re.search(r'\bI\s*=\s*(\d+)', line) npoints = int(m_i.group(1)) if m_i else None i += 1 xc_list, cp_list = [], [] count = 0 while i < len(lines) and (npoints is None or count < npoints): data_line = lines[i].strip() if not data_line: i += 1 continue if data_line.upper().startswith(("ZONE", "TITLE", "VARIABLES")): break parts = data_line.split() try: vals = [float(p) for p in parts] xc_list.append(vals[2]) # column 2: X/L (= X/C) cp_list.append(vals[4]) # column 4: CP count += 1 except (ValueError, IndexError): pass i += 1 zones.append({ "label": label, "xc": np.array(xc_list), "cp": np.array(cp_list), }) else: i += 1 return zones # --------------------------------------------------------------------------- # Plotting # --------------------------------------------------------------------------- def plot_sections(exp_zones, lava_zones, output_dir: Path): output_dir.mkdir(parents=True, exist_ok=True) n_sections = min(len(exp_zones), len(lava_zones), len(YB_VALUES)) for idx in range(n_sections): yb = YB_VALUES[idx] exp = exp_zones[idx] lava = lava_zones[idx] fig, ax = plt.subplots(figsize=(7, 5)) # Experiment: black filled squares ax.plot(exp["xc"], exp["cp"], marker="s", markersize=6, linestyle="none", color="black", markerfacecolor="black", label="Experiment") # LAVA: red line — plot in original data order so upper/lower # surface points are not interleaved when connecting the line ax.plot(lava["xc"], lava["cp"], color="red", linewidth=1.5, label="LAVA (unst)") # Invert y-axis so negative Cp is at the top ax.invert_yaxis() ax.set_xlabel("x/c", fontsize=12) ax.set_ylabel(r"$C_p$", fontsize=12) ax.set_title(rf"ONERA M6 — $y/b = {yb}$", fontsize=13) ax.legend(fontsize=11) ax.grid(True, linestyle="--", alpha=0.4) fname = output_dir / f"Cp_section_{idx+1}_yb{str(yb).replace('.','p')}.png" fig.savefig(fname, dpi=150, bbox_inches="tight") plt.close(fig) print(f" Saved: {fname}") # --------------------------------------------------------------------------- # Main # --------------------------------------------------------------------------- def main(): here = Path(__file__).parent parser = argparse.ArgumentParser(description="Plot Cp sections for ONERA M6") parser.add_argument("--exp", default=str(here / "exp_case_2308.dat"), help="Path to experimental data file") parser.add_argument("--lava", default=str(here / "surface" / "Unst.cutdata.000000500.dat"), help="Path to LAVA cutdata file") parser.add_argument("--out", default=str(here / "cp_plots"), help="Output directory for PNG figures") args = parser.parse_args() exp_file = Path(args.exp) lava_file = Path(args.lava) out_dir = Path(args.out) print(f"Experiment : {exp_file}") print(f"LAVA data : {lava_file}") print(f"Output dir : {out_dir}") exp_zones = parse_experiment(exp_file) lava_zones = parse_lava_cutdata(lava_file) print(f"\nParsed {len(exp_zones)} experimental zones") print(f"Parsed {len(lava_zones)} LAVA zones") plot_sections(exp_zones, lava_zones, out_dir) print("\nDone.") if __name__ == "__main__": main()