SheCodes Lab

# Find out where you are (your current folder)
pwd              # Mac/Linux
cd               # Windows (just typing cd shows your location)

# List files in the current folder
ls               # Mac/Linux
dir              # Windows

# Move into a folder
cd Documents
cd Documents/my_code

# Go back up one folder
cd ..

# Run a Python file
python hello.py
python3 hello.py  # if the above doesn't work

# Install a Python library
pip install numpy
pip install pandas matplotlib scipy

# 1. Open VS Code
# 2. File → New File → save it as:  hello.py
#    (the .py tells the computer it's Python)

# 3. Type this into the file:
print("Hello! I am learning to code.")
print("This is my first program.")

# 4. Open the terminal in VS Code: View → Terminal
# 5. Type:  python hello.py  and press Enter
# 6. You should see:
#    Hello! I am learning to code.
#    This is my first program.

# Congratulations — you just ran your first program.

import numpy  as np
import pandas as pd
import matplotlib.pyplot as plt

# 1. Load test data from the testing machine CSV export
df = pd.read_csv("al6061_tensile.csv")           # columns: strain_pct, stress_MPa

# 2. Find the 0.2% proof stress (yield strength)
offset_line = 0.2 + df["strain_pct"] * 0         # horizontal offset
elastic     = 68900 * (df["strain_pct"] - 0.2) / 100   # E = 68.9 GPa

# 3. Compute basic stats
uts = df["stress_MPa"].max()
print(f"UTS:              {uts:.1f} MPa")
print(f"Fracture strain:  {df['strain_pct'].iloc[-1]:.2f}%")

# 4. Plot the stress-strain curve
plt.figure(figsize=(9, 6))
plt.plot(df["strain_pct"], df["stress_MPa"],
         color="steelblue", linewidth=2, label="Al 6061-T6")
plt.axhline(uts, color="red", linestyle="--", label=f"UTS = {uts:.0f} MPa")
plt.xlabel("Strain (%)");  plt.ylabel("Stress (MPa)")
plt.title("Tensile Test — Al 6061-T6")
plt.legend();  plt.grid(True, alpha=0.3)
plt.savefig("stress_strain.png", dpi=300)
plt.show()

import numpy as np

# Create arrays from data
stress  = np.array([0, 50, 100, 150, 200, 250])   # MPa
strain  = np.array([0, 0.07, 0.14, 0.22, 0.29, 0.36]) # %

# Maths on every element at once — no loop
stress_psi = stress * 145.038    # convert MPa → psi (all 6 values at once)
strain_dec = strain / 100         # percent → decimal

# Generate engineering sequences
angles    = np.linspace(0, 360, 361)   # 0° to 360°, every 1°
time      = np.linspace(0, 10, 1000)   # 0 to 10 seconds, 1000 points
thicknesses = np.arange(1, 20, 0.5)    # 1mm to 19.5mm, step 0.5mm

# Array properties
print(stress.shape)    # (6,) — 1D array with 6 elements
print(stress.size)     # 6
print(stress.dtype)    # int64

import numpy as np

# Trig — angles in radians (convert first)
aoa_deg = np.array([0, 2, 4, 6, 8, 10, 12])    # angle of attack, degrees
aoa_rad = np.radians(aoa_deg)

# Thin aerofoil theory: CL ≈ 2π·sin(α)  (linearised)
CL = 2 * np.pi * np.sin(aoa_rad)
print("CL values:", np.round(CL, 3))

# Decompose a velocity vector into components
V_total = 250       # m/s (TAS)
gamma   = 15        # flight path angle, degrees
Vx = V_total * np.cos(np.radians(gamma))   # horizontal component
Vz = V_total * np.sin(np.radians(gamma))   # vertical component (climb)
print(f"Vx = {Vx:.1f} m/s,  Vz = {Vz:.1f} m/s")

# 2D rotation matrix — rotate a force vector by θ
theta = np.radians(30)
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])
F = np.array([1000, 0])       # 1000 N in x direction
F_rotated = R @ F               # @ is matrix multiply
print(f"Rotated force: {F_rotated.round(1)} N")

import numpy as np

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
# Example: 3-bar truss — solve for nodal displacements
# K·u = F  →  stiffness matrix × displacements = forces
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

K = np.array([[ 3, -1,  0],
              [-1,  2, -1],
              [ 0, -1,  1]], dtype=float) # stiffness matrix

F = np.array([0, 10000, 0])               # applied forces, N

# Solve — equivalent to MATLAB's u = K\F
u = np.linalg.solve(K, F)
print("Displacements (mm):", u.round(4))

# Matrix properties engineers need
print(f"Determinant:   {np.linalg.det(K):.2f}")
print(f"Condition no.: {np.linalg.cond(K):.2f}")  # ill-conditioned if large

eigenvalues, eigenvectors = np.linalg.eig(K)
print("Natural frequencies (proportional):", np.sqrt(eigenvalues).round(3))

import numpy as np

# ISA constants
T0   = 288.15   # sea level temperature, K
P0   = 101325   # sea level pressure, Pa
rho0 = 1.225    # sea level density, kg/m³
L    = 0.0065   # lapse rate, K/m
g    = 9.80665  # gravitational acceleration
R    = 287.05   # specific gas constant, J/(kg·K)

# Altitude range: sea level to 11 km (troposphere), 100 points
h = np.linspace(0, 11000, 100)     # metres

# ISA equations — applied to entire array at once
T   = T0 - L * h                         # temperature, K
P   = P0 * (T / T0) ** (g / (L * R))    # pressure, Pa
rho = P / (R * T)                         # density, kg/m³
a   = np.sqrt(1.4 * R * T)               # speed of sound, m/s

# Print at specific altitudes
for alt in [0, 3000, 6000, 9000, 11000]:
    idx = np.argmin(np.abs(h - alt))     # find nearest index
    print(f"h={alt:5d}m  T={T[idx]:.1f}K  rho={rho[idx]:.4f}kg/m³  a={a[idx]:.1f}m/s")

import numpy as np

# UTS results from 20 tensile specimens of Ti-6Al-4V (MPa)
uts = np.array([940, 952, 938, 961, 945, 958, 933, 970,
                948, 955, 942, 964, 950, 937, 956, 947,
                962, 943, 959, 951])

mean   = np.mean(uts)
std    = np.std(uts, ddof=1)        # ddof=1 for sample std dev
cv     = (std / mean) * 100         # coefficient of variation, %
se     = std / np.sqrt(len(uts))    # standard error

print(f"Mean UTS:           {mean:.1f} MPa")
print(f"Std deviation:      {std:.1f} MPa")
print(f"Coefficient of var: {cv:.2f}%")
print(f"Min / Max:          {uts.min()} / {uts.max()} MPa")
print(f"95% CI (approx):    {mean:.1f} ± {2*se:.1f} MPa")

# Spot outliers: values beyond 2 standard deviations
outliers = uts[np.abs(uts - mean) > 2 * std]
print(f"Outliers:           {outliers}")

import numpy as np
import matplotlib.pyplot as plt

# Fatigue data: stress amplitude (MPa) vs log10(cycles to failure)
stress   = np.array([400, 350, 300, 280, 260, 240, 220, 200])
log_N    = np.array([3.8, 4.2, 4.9, 5.3, 5.7, 6.1, 6.6, 7.0])

# Linear regression: fit a straight line y = m*x + c
coeffs  = np.polyfit(log_N, stress, 1)   # degree 1 = linear
m, c    = coeffs
print(f"Slope (m):     {m:.2f}")
print(f"Intercept (c): {c:.2f}")
print(f"Equation:      stress = {m:.1f} * log10(N) + {c:.1f}")

# Calculate R² (coefficient of determination)
predicted  = np.polyval(coeffs, log_N)
ss_res     = np.sum((stress - predicted) ** 2)
ss_tot     = np.sum((stress - stress.mean()) ** 2)
r_squared  = 1 - (ss_res / ss_tot)
print(f"R² = {r_squared:.4f}")           # closer to 1.0 = better fit

# Predict: what stress gives 10^6 cycles?
stress_at_1M = np.polyval(coeffs, 6.0)
print(f"Fatigue limit (10^6 cycles): {stress_at_1M:.1f} MPa")

# Plot
fit_line = np.linspace(3.5, 7.5, 100)
plt.scatter(log_N, stress, color="steelblue", s=60, label="Test data", zorder=5)
plt.plot(fit_line, np.polyval(coeffs, fit_line), "r--",
         label=f"Linear fit  R²={r_squared:.3f}")
plt.xlabel("log₁₀(N)");  plt.ylabel("Stress amplitude (MPa)")
plt.title("S-N Curve — Al 2024-T3")
plt.legend();  plt.grid(True, alpha=0.3)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Experimental lift coefficient data (clean wing, subsonic)
alpha = np.array([-4, -2,  0,  2,  4,  6,  8, 10, 12])  # AoA, deg
CL    = np.array([-0.28, -0.06, 0.18, 0.42, 0.65, 0.88,
                    1.05, 1.18, 1.22])

# Fit a 2nd-degree polynomial (captures the slight nonlinearity)
coeffs  = np.polyfit(alpha, CL, 2)
print(f"CL = {coeffs[0]:.5f}·α² + {coeffs[1]:.4f}·α + {coeffs[2]:.4f}")

# Lift slope dCL/dα at α=0° (per degree)
lift_slope = 2 * coeffs[0] * 0 + coeffs[1]
print(f"Lift slope: {lift_slope:.4f} /deg  ({lift_slope * 57.3:.3f} /rad)")

# Zero-lift angle
zero_lift_alpha = np.roots(coeffs)[1]  # find root of polynomial
print(f"Zero-lift angle: {zero_lift_alpha:.2f}°")

# Plot the fit against the data
alpha_smooth = np.linspace(-5, 14, 200)
CL_fit = np.polyval(coeffs, alpha_smooth)
plt.scatter(alpha, CL, color="steelblue", s=60, zorder=5, label="Wind tunnel data")
plt.plot(alpha_smooth, CL_fit, "r-", linewidth=2, label="Quadratic fit")
plt.axhline(0, color="gray", lw=0.5)
plt.axvline(0, color="gray", lw=0.5)
plt.xlabel("Angle of Attack α (°)")
plt.ylabel("Lift Coefficient CL")
plt.title("Lift Polar — Clean Wing")
plt.legend();  plt.grid(True, alpha=0.3)
plt.show()

import pandas as pd

# Load CSV exported from test rig or DAS (data acquisition system)
df = pd.read_csv("wing_fatigue_tests.csv")

# Or load from Excel directly
df = pd.read_excel("material_test_results.xlsx", sheet_name="Fatigue")

# First look — always do this before any analysis
print(df.head(5))         # first 5 rows
print(df.shape)           # (rows, columns)
print(df.columns.tolist()) # column names
print(df.dtypes)          # data types — check nothing read as text
print(df.isnull().sum())  # missing values per column
print(df.describe())      # count, mean, std, min, 25%, 50%, 75%, max

# Select a single column
cycles = df["cycles_to_failure"]

# Filter rows: only specimens that failed (not run-outs)
failed  = df[df["status"] == "FAILURE"]
runouts = df[df["status"] == "RUNOUT"]
print(f"Failures: {len(failed)},  Run-outs: {len(runouts)}")

# Filter: specific material AND stress range
al_high = df[
    (df["material"] == "Al2024-T3") &
    (df["stress_MPa"] >= 200) &
    (df["stress_MPa"] <= 350)
]

# Filter to several materials at once
metals = df[df["material"].isin(["Al2024-T3", "Ti6Al4V", "Steel4340"])]

# Sort by cycles descending
df_sorted = df.sort_values("cycles_to_failure", ascending=False)

# Compute stress from raw load and cross-section area
df["stress_MPa"]   = df["load_N"]  / df["area_mm2"]

# Compute engineering strain from displacement and gauge length
df["strain_pct"]   = (df["disp_mm"] / df["gauge_mm"]) * 100

# Flag specimens that exceeded design limit
df["overload"]     = df["stress_MPa"] > 250

# Mean stress and stress ratio for fatigue
df["stress_mean"]  = (df["stress_max"] + df["stress_min"]) / 2
df["stress_ratio"] = df["stress_min"] / df["stress_max"]

# Group by material — get mean, std, count per group (like a pivot table)
summary = df.groupby("material").agg(
    mean_UTS   = ("stress_MPa", "mean"),
    std_UTS    = ("stress_MPa", "std"),
    n_specimens= ("stress_MPa", "count")
).round(2)
print(summary)

import pandas as pd
import numpy  as np

# 1. Load raw test data
df = pd.read_csv("fatigue_database.csv")
print(f"Loaded {len(df)} records.")

# 2. Clean — remove rows with missing critical values
df = df.dropna(subset=["stress_MPa", "cycles_to_failure"])

# 3. Compute log cycles (for S-N plot)
df["log_N"] = np.log10(df["cycles_to_failure"])

# 4. Filter to material of interest
mat = df[df["material"] == "Al2024-T3"]

# 5. Linear regression on S-N data
coeffs = np.polyfit(mat["log_N"], mat["stress_MPa"], 1)
print(f"Fatigue limit (10^7): {np.polyval(coeffs, 7):.1f} MPa")

# 6. Export processed results
mat.to_csv("al2024_processed.csv", index=False)
print("Saved al2024_processed.csv")

import matplotlib.pyplot as plt
import numpy as np

# Simulated tensile test data — Al 6061-T6
strain = np.array([0,    0.1,  0.2,  0.3,  0.45, 0.7,
                    1.0,  1.5,  2.0,  3.0,  4.5,  6.0,  8.0])
stress = np.array([0,    69,  138,  207,  270, 290,
                    300,  310,  318,  325,  310,  285,  240])

uts_idx   = np.argmax(stress)
uts_val   = stress[uts_idx]
uts_strain= strain[uts_idx]

fig, ax = plt.subplots(figsize=(9, 6))

# Main curve
ax.plot(strain, stress, color="#2563EB", linewidth=2.5, label="Al 6061-T6")

# Elastic modulus line (first two points)
ax.plot([0, 0.3], [0, 68.9*0.003*1000], "k--", lw=1, label="Elastic region")

# Annotate UTS
ax.annotate(f"UTS = {uts_val} MPa",
            xy=(uts_strain, uts_val),
            xytext=(uts_strain - 1.5, uts_val - 25),
            arrowprops=dict(arrowstyle="->", color="red"),
            color="red", fontsize=10)

ax.set_xlabel("Engineering Strain (%)", fontsize=12)
ax.set_ylabel("Engineering Stress (MPa)", fontsize=12)
ax.set_title("Tensile Test — Al 6061-T6", fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("stress_strain.png", dpi=300)
plt.show()

import matplotlib.pyplot as plt
import numpy as np

t    = np.linspace(0, 120, 1200)          # 0–120 s at 10 Hz
alt  = 1000 * (t/120)**0.5                # altitude climb
spd  = 60 + 80 * (t/120)                 # speed ramp
aoa  = 8 * np.exp(-t/30) + 2              # AoA reducing after rotation
nz   = 1 + 0.3 * np.sin(2*np.pi*t/15)    # normal load factor

fig, axes = plt.subplots(4, 1, figsize=(11, 10), sharex=True)
fig.suptitle("Flight Test Data — Climb Segment", fontsize=14, y=0.98)

data   = [alt, spd, aoa, nz]
labels = ["Altitude (m)", "Airspeed (m/s)", "AoA (°)", "Load Factor Nz"]
colors = ["#2563EB", "#16A34A", "#DC2626", "#9333EA"]

for ax, d, lbl, col in zip(axes, data, labels, colors):
    ax.plot(t, d, color=col, linewidth=1.5)
    ax.set_ylabel(lbl, fontsize=10)
    ax.grid(True, alpha=0.25)

axes[3].axhline(2.5, color="red", ls="--", lw=1, label="Limit load")
axes[3].legend(fontsize=9)
axes[3].set_xlabel("Time (s)", fontsize=11)

plt.tight_layout()
plt.savefig("flight_data.png", dpi=300)
plt.show()

import matplotlib.pyplot as plt
import numpy as np

# Simulated dataset: fatigue tests at different temperatures
np.random.seed(42)
n       = 80
stress  = np.random.uniform(150, 400, n)
log_N   = 8.5 - 0.012 * stress + np.random.normal(0, 0.2, n)
temp    = np.random.uniform(20, 300, n)    # temperature, °C

fig, ax = plt.subplots(figsize=(9, 6))

sc = ax.scatter(log_N, stress, c=temp, cmap="plasma", s=50, alpha=0.8)
cbar = plt.colorbar(sc, ax=ax)
cbar.set_label("Temperature (°C)", fontsize=11)

ax.set_xlabel("log₁₀(N) — Cycles to Failure", fontsize=12)
ax.set_ylabel("Stress Amplitude (MPa)",       fontsize=12)
ax.set_title("S-N Data Coloured by Test Temperature", fontsize=13)
ax.grid(True, alpha=0.25)
plt.tight_layout()
plt.savefig("sn_temperature.png", dpi=300)
plt.show()

import pandas as pd

# Read first sheet (default)
df = pd.read_excel("test_results.xlsx")

# Read a named sheet
df = pd.read_excel("test_results.xlsx", sheet_name="Fatigue")

# Read all sheets at once (returns a dict)
all_sheets = pd.read_excel("test_results.xlsx", sheet_name=None)
for name, sheet_df in all_sheets.items():
    print(f"Sheet '{name}': {len(sheet_df)} rows")

# Skip metadata rows at the top (common in lab exports)
df = pd.read_excel("test_results.xlsx", skiprows=4, header=0)

# Read only specific columns to save memory
df = pd.read_excel("test_results.xlsx",
                   usecols=["specimen_id", "material", "UTS_MPa", "cycles"])

import pandas as pd
import numpy  as np

# Your analysed data
raw_df = pd.read_csv("fatigue_tests.csv")
raw_df["log_N"] = np.log10(raw_df["cycles_to_failure"])

# Summary statistics per material
summary = raw_df.groupby("material").agg(
    n_specimens = ("UTS_MPa", "count"),
    mean_UTS    = ("UTS_MPa", "mean"),
    std_UTS     = ("UTS_MPa", "std"),
    min_cycles  = ("cycles_to_failure", "min"),
    max_cycles  = ("cycles_to_failure", "max")
).round(2)

# Write to Excel with multiple sheets
with pd.ExcelWriter("fatigue_report.xlsx", engine="openpyxl") as writer:
    raw_df.to_excel(writer, sheet_name="Raw Data",   index=False)
    summary.to_excel(writer, sheet_name="Summary",    index=True)

print("Saved fatigue_report.xlsx with 2 sheets.")

import pandas as pd
import numpy  as np
import glob
import os

# Find every .xlsx file in the test_data folder
files = glob.glob("test_data/*.xlsx")
print(f"Found {len(files)} test files.")

results = []

for filepath in files:
    specimen_id = os.path.basename(filepath).replace(".xlsx", "")

    # Load each file
    df = pd.read_excel(filepath)

    # Compute stress from raw columns
    df["stress_MPa"] = df["load_N"] / df["area_mm2"]

    # Extract summary for this specimen
    row = {
        "specimen":    specimen_id,
        "material":    df["material"].iloc[0],
        "UTS_MPa":     df["stress_MPa"].max().round(1),
        "E_GPa":       (df["stress_MPa"].iloc[5] / df["strain_pct"].iloc[5] * 100).round(1),
        "elong_pct":   df["strain_pct"].iloc[-1]
    }
    results.append(row)
    print(f"  Processed {specimen_id}  →  UTS = {row['UTS_MPa']} MPa")

# Compile all results
report = pd.DataFrame(results)
print("\n=== BATCH SUMMARY ===")
print(report)

# Append overall stats to the bottom
report.to_excel("batch_summary.xlsx", index=False)
print("\nSaved batch_summary.xlsx")

import numpy as np
from scipy import fft
import matplotlib.pyplot as plt

# Simulate accelerometer data: 50 Hz + 120 Hz signal + noise
fs  = 1000                              # sample rate, Hz
t   = np.linspace(0, 1, fs, endpoint=False) # 1 second
sig = (np.sin(2*np.pi*50*t)          # 50 Hz component
     + 0.5*np.sin(2*np.pi*120*t)      # 120 Hz component
     + 0.3*np.random.randn(fs))        # noise

# Compute FFT
N      = len(sig)
yf     = fft.fft(sig)
xf     = fft.fftfreq(N, 1/fs)          # frequency axis
power  = 2/N * np.abs(yf[:N//2])       # one-sided power spectrum
freqs  = xf[:N//2]

# Find dominant frequency
peak_freq = freqs[np.argmax(power)]
print(f"Dominant frequency: {peak_freq:.1f} Hz")

# Plot spectrum
plt.figure(figsize=(10, 4))
plt.plot(freqs, power, color="steelblue", lw=1.2)
plt.xlabel("Frequency (Hz)");  plt.ylabel("Amplitude")
plt.title("FFT Power Spectrum");  plt.grid(True, alpha=0.3)
plt.xlim(0, 300)
plt.show()

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Define the equation to fit — power law: y = a * x^b
def power_law(x, a, b):
    return a * x ** b

# Creep strain data: time (hours) vs strain (%)
time   = np.array([1, 5, 10, 50, 100, 500, 1000])
strain = np.array([0.12, 0.21, 0.28, 0.48, 0.61, 1.05, 1.32])

# Fit — popt = best parameters, pcov = covariance (uncertainty)
popt, pcov = curve_fit(power_law, time, strain)
a, b = popt
perr = np.sqrt(np.diag(pcov))  # 1-sigma uncertainty

print(f"a = {a:.4f} ± {perr[0]:.4f}")
print(f"b = {b:.4f} ± {perr[1]:.4f}")
print(f"Equation: strain = {a:.4f} * t^{b:.4f}")

# Predict strain at 2000 hours
print(f"Predicted strain at 2000h: {power_law(2000, *popt):.3f}%")

# Plot
t_smooth = np.logspace(0, 4, 200)
plt.scatter(time, strain, s=60, zorder=5, label="Data")
plt.plot(t_smooth, power_law(t_smooth, *popt), "r-", label=f"Fit: {a:.3f}·t^{b:.3f}")
plt.xscale("log");  plt.yscale("log")
plt.xlabel("Time (h)");  plt.ylabel("Strain (%)")
plt.title("Creep Curve — Power Law Fit")
plt.legend();  plt.grid(True, alpha=0.3)
plt.show()

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
g0    = 9.80665      # m/s²
m0    = 10000        # initial mass, kg
mdot  = 50           # mass flow rate, kg/s
Isp   = 300          # specific impulse, s
ve    = Isp * g0     # exhaust velocity, m/s
Cd    = 0.3          # drag coefficient
A     = 1.0          # reference area, m²
rho0  = 1.225        # sea level density, kg/m³

def rocket_eom(t, y):
    """State vector: y = [altitude, velocity, mass]"""
    h, v, m = y
    rho  = rho0 * np.exp(-h / 8500)       # exponential atmosphere
    drag = 0.5 * rho * v**2 * Cd * A
    thrust = mdot * ve if m > m0 - mdot*60 else 0   # burn for 60 s
    dhdt = v
    dvdt = (thrust - drag) / m - g0
    dmdt = -mdot if thrust > 0 else 0
    return [dhdt, dvdt, dmdt]

# Solve from t=0 to t=120s
sol = solve_ivp(rocket_eom, [0, 120], [0, 0, m0],
                max_step=0.5, dense_output=True)

t   = sol.t
h   = sol.y[0] / 1000    # altitude, km
v   = sol.y[1]             # velocity, m/s
print(f"Max altitude: {h.max():.2f} km")
print(f"Max velocity: {v.max():.1f} m/s")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
ax1.plot(t, h, color="#2563EB");  ax1.set_ylabel("Altitude (km)")
ax2.plot(t, v, color="#DC2626");  ax2.set_ylabel("Velocity (m/s)")
for ax in (ax1, ax2):
    ax.set_xlabel("Time (s)");  ax.grid(True, alpha=0.3)
plt.suptitle("Rocket Trajectory");  plt.tight_layout()
plt.show()

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

fs     = 500                          # sample rate, Hz
t      = np.linspace(0, 2, 2*fs)     # 2 seconds
clean  = np.sin(2*np.pi*5*t)          # 5 Hz true signal
noisy  = clean + 0.5*np.random.randn(len(t))

# Design a 4th-order Butterworth low-pass filter, cutoff = 20 Hz
b, a   = signal.butter(4, 20, btype="low", fs=fs)
filtered = signal.filtfilt(b, a, noisy)    # zero-phase filtering

plt.figure(figsize=(10, 4))
plt.plot(t, noisy,    alpha=0.4, label="Noisy",    lw=1)
plt.plot(t, filtered, color="#DC2626", label="Filtered", lw=2)
plt.xlabel("Time (s)");  plt.ylabel("Signal")
plt.title("Butterworth Low-Pass Filter")
plt.legend();  plt.grid(True, alpha=0.3)
plt.show()

import numpy as np

# BUG 1: Integer division gives wrong answer silently
E   = 200000   # MPa
eps = 3        # strain — but these are both ints!
stress = E * eps / 1000    # Python 3: gives 0.6 ✓
# In Python 2 or C++: 200000 * 3 / 1000 = 600, not 0.6
# In C++: int/int = int — use 1000.0 to force float division

# BUG 2: Comparing floats with == (almost always wrong)
a = 0.1 + 0.2
print(a == 0.3)              # False! (floating point)
print(np.isclose(a, 0.3))   # True  ← use this instead

# BUG 3: Mutating a list you think you're copying
original = [1, 2, 3]
copy     = original          # NOT a copy — same object!
copy[0]  = 99
print(original)              # [99, 2, 3] — original changed!
# Fix:
copy     = original.copy()  # or list(original) or original[:]
# NumPy:
arr2 = arr.copy()           # NOT arr2 = arr

# BUG 4: Degrees vs radians (kills calculations silently)
angle = 45                     # degrees
print(np.sin(angle))          # 0.8509 — WRONG, treats as radians
print(np.sin(np.radians(45))) # 0.7071 — correct

# BUG 5: NaN propagation — one bad value corrupts everything
data = np.array([1.0, 2.0, np.nan, 4.0])
print(np.mean(data))     # nan — mean fails silently
print(np.nanmean(data))  # 2.333 — use nan-safe versions
print(np.any(np.isnan(data)))  # True — always check your data

import time
import numpy as np

# Simple timing with perf_counter (high resolution)
start = time.perf_counter()

# --- your code here ---
result = np.sum(np.random.randn(1_000_000))

elapsed = time.perf_counter() - start
print(f"Elapsed: {elapsed*1000:.2f} ms")

# --- Python loop vs NumPy --- why it matters
data = np.random.randn(100_000)

# Slow: Python loop
t0 = time.perf_counter()
total = 0
for x in data:
    total += x ** 2
print(f"Loop:  {(time.perf_counter()-t0)*1000:.1f} ms")

# Fast: NumPy vectorised
t0 = time.perf_counter()
total = np.sum(data ** 2)
print(f"NumPy: {(time.perf_counter()-t0)*1000:.1f} ms")
# NumPy is typically 50-200x faster on large arrays

# Full profiling with cProfile (run from terminal)
# python -m cProfile -s cumulative your_script.py

import numpy as np

def compute_stress(force_N, area_mm2):
    # Guard against obviously wrong inputs
    assert force_N >= 0,    "Force must be non-negative"
    assert area_mm2 > 0,   "Area must be positive"
    stress = force_N / area_mm2
    assert stress < 2000, f"Stress {stress:.1f} MPa exceeds material limit"
    return stress

# Validate array shapes before matrix operations
A = np.ones((3, 3))
b = np.ones(3)
assert A.shape[1] == b.shape[0], "Matrix dimensions must match"

# Check no NaN crept in after operations
result = np.linalg.solve(A, b)
assert not np.any(np.isnan(result)), "NaN in result — check input matrix"
print("All checks passed.")

from scipy import constants as c
import numpy as np

# Physical constants (exact scipy values)
print(c.g)          # 9.80665    m/s²  standard gravity
print(c.R)          # 8.314...   J/(mol·K) universal gas constant
print(c.atm)        # 101325     Pa    standard atmosphere
print(c.speed_of_light) # 299792458 m/s
print(c.pi)         # same as np.pi

# Aerospace / mechanical constants — define these at top of every script
g0    = 9.80665     # m/s²  — standard gravity (ISA)
R_air = 287.058    # J/(kg·K) — specific gas constant for dry air
gamma = 1.4        # —  ratio of specific heats for air
P_sl  = 101325.0   # Pa   — sea level pressure (ISA)
T_sl  = 288.15     # K    — sea level temperature (ISA)
rho_sl= 1.225      # kg/m³ — sea level density (ISA)
a_sl  = np.sqrt(gamma * R_air * T_sl) # 340.3 m/s — speed of sound

# Unit conversions — name them clearly
MPa_to_psi  = 145.038
psi_to_MPa  = 1 / MPa_to_psi
kts_to_ms   = 0.514444
ms_to_kts   = 1 / kts_to_ms
ft_to_m     = 0.3048
m_to_ft     = 1 / ft_to_m
lbf_to_N    = 4.44822
deg_to_rad  = np.pi / 180    # or use np.radians()

# Fixed decimal places
stress = 247.38291
print(f"{stress:.2f}")        # 247.38        (2 decimal places)
print(f"{stress:.0f}")        # 247            (round to integer)

# Significant figures
val = 0.0001234567
print(f"{val:.3g}")           # 0.000123       (3 sig figs)
print(f"{val:.4g}")           # 0.0001235

# Scientific notation
Re = 4_500_000
print(f"{Re:.2e}")            # 4.50e+06
print(f"{Re:,.0f}")           # 4,500,000      (thousands separator)

# Width padding — align columns
data = [("Steel", 400, 210), ("Aluminium", 270, 69), ("CFRP", 600, 150)]
print(f"{'Material':12} {'UTS (MPa)':>10} {'E (GPa)':>10}")
print("-" * 34)
for name, uts, E in data:
    print(f"{name:12} {uts:>10.0f} {E:>10.0f}")

from pint import UnitRegistry

ureg = UnitRegistry()
Q    = ureg.Quantity

# Attach units to values
force  = Q(10000, "N")
area   = Q(50,    "mm**2")

# Arithmetic is unit-aware
stress = force / area
print(stress)               # 200.0 N/mm²
print(stress.to("MPa"))     # 200.0 MPa  (N/mm² = MPa)
print(stress.to("psi"))     # 29007.5 psi

# Unit mismatch raises error immediately
height = Q(100, "ft")
# force + height  →  DimensionalityError: cannot add [force] and [length]

# Convert velocity
speed = Q(250, "m/s")
print(speed.to("knots"))   # 486.0 knots
print(speed.to("mph"))     # 559.2 mph

import numpy as np

# Cl vs AoA lookup table (from wind tunnel data)
aoa_data = np.array([ 0,    2,    4,    6,    8,   10,   12])
CL_data  = np.array([ 0.18, 0.42, 0.65, 0.88, 1.05, 1.18, 1.22])

# Interpolate at arbitrary angle
aoa_query = 5.5
CL_at_5p5 = np.interp(aoa_query, aoa_data, CL_data)
print(f"CL at {aoa_query}°: {CL_at_5p5:.4f}")

# Interpolate multiple angles at once
query_angles = np.array([1, 3, 5, 7, 9, 11])
CL_interp    = np.interp(query_angles, aoa_data, CL_data)
print(CL_interp)

# scipy for more control (cubic spline — smoother)
from scipy.interpolate import CubicSpline
cs = CubicSpline(aoa_data, CL_data)
print(cs(5.5))                   # smoother than linear interp
print(cs(5.5, 1))               # derivative: dCL/dα at 5.5°

import numpy as np

# ISA constants
T0    = 288.15     # K    sea-level temperature
P0    = 101325.0   # Pa   sea-level pressure
rho0  = 1.225      # kg/m³ sea-level density
L     = 0.0065     # K/m  troposphere lapse rate
g0    = 9.80665    # m/s²
R     = 287.058    # J/(kg·K) dry air
gamma = 1.4
H_trop = 11000     # m — tropopause altitude
T_trop = T0 - L * H_trop   # 216.65 K at tropopause
P_trop = P0 * (T_trop / T0) ** (g0 / (L * R))

def isa(h_m):
    """
    ISA atmosphere for scalar or array altitude h_m (metres).
    Returns: T (K), P (Pa), rho (kg/m³), a (m/s)
    """
    h = np.atleast_1d(np.asarray(h_m, dtype=float))
    T   = np.where(h <= H_trop,
                  T0 - L * h,
                  T_trop)
    P   = np.where(h <= H_trop,
                  P0 * (T / T0) ** (g0 / (L * R)),
                  P_trop * np.exp(-g0 * (h - H_trop) / (R * T_trop)))
    rho = P / (R * T)
    a   = np.sqrt(gamma * R * T)
    return T.squeeze(), P.squeeze(), rho.squeeze(), a.squeeze()

# --- Usage ---
altitudes_ft = np.array([0, 10000, 20000, 30000, 35000, 39000])
altitudes_m  = altitudes_ft * 0.3048

T, P, rho, a = isa(altitudes_m)

print(f"{'Alt (ft)':10} {'T (K)':8} {'P (hPa)':10} {'rho (kg/m³)':13} {'a (kts)'}")
print("-"*55)
for i, alt in enumerate(altitudes_ft):
    print(f"{alt:10.0f} {T[i]:8.2f} {P[i]/100:10.2f} {rho[i]:13.4f} {a[i]/0.514444:7.1f}")

def mach_to_tas(M, h_m):
    """True Airspeed from Mach and altitude"""
    T, _, _, a = isa(h_m)
    return M * a                             # m/s

def tas_to_eas(tas, h_m):
    """Equivalent Airspeed = TAS * sqrt(sigma)"""
    _, _, rho, _ = isa(h_m)
    sigma = rho / rho0
    return tas * np.sqrt(sigma)

def cas_to_tas(cas, h_m):
    """CAS to TAS using compressibility correction"""
    T, P, rho, _ = isa(h_m)
    # Subsonic isentropic relation
    qc = P0 * ((1 + 0.2 * (cas / (isa(0)[3]))**2) ** 3.5 - 1)
    M   = np.sqrt(5 * ((qc / P + 1) ** (2/7) - 1))
    _, _, _, a = isa(h_m)
    return M * a

# Example: A320 cruise at FL350, M0.78
h_cruise = 35000 * 0.3048      # FL350 in metres
M        = 0.78
TAS      = mach_to_tas(M, h_cruise)
EAS      = tas_to_eas(TAS, h_cruise)

print(f"FL350, M{M}:")
print(f"  TAS = {TAS:.1f} m/s  ({TAS/0.514444:.1f} kts)")
print(f"  EAS = {EAS:.1f} m/s  ({EAS/0.514444:.1f} kts)")

import numpy as np
import matplotlib.pyplot as plt

# Aircraft parameters (A320-like)
W       = 65000    # kg  cruise weight
S       = 122.6    # m²  wing area
CD0     = 0.0240   # zero-lift drag
k       = 0.0375   # induced drag factor (1/π·AR·e)
eta     = 0.30     # overall propulsive efficiency (cruise)
LHV     = 43.2e6   # J/kg  Jet-A lower heating value

# Cruise conditions
h_m     = 35000 * 0.3048
T, P, rho, a = isa(h_m)

# Mach sweep
M_range = np.linspace(0.65, 0.84, 200)
V       = M_range * a                             # TAS, m/s
q       = 0.5 * rho * V**2                        # dynamic pressure
CL      = W * 9.80665 / (q * S)                  # level flight CL
CD      = CD0 + k * CL**2                         # total drag
D       = q * S * CD                               # drag force, N
ff_kgs  = D * V / (eta * LHV)                    # fuel flow, kg/s
ff_kghr = ff_kgs * 3600                             # fuel flow, kg/hr
SR      = V / ff_kgs / 1000                       # specific range, km/kg

# Cost Index analysis
# CI units: kg/min  (fuel cost equivalent of 1 min of time)
def opt_mach_for_ci(CI_kgmin):
    """Find Mach that minimises total cost for given CI"""
    time_cost_per_km = CI_kgmin / (V * 60 / 1000)   # kg_eq/km from time
    fuel_cost_per_km = ff_kgs / (V / 1000)           # kg/km
    total_cost       = time_cost_per_km + fuel_cost_per_km
    idx              = np.argmin(total_cost)
    return M_range[idx], ff_kghr[idx]

print(f"{'CI (kg/min)':14} {'Opt Mach':12} {'Fuel Flow (kg/hr)'}")
print("-"*42)
for CI in [0, 10, 20, 40, 60, 80, 100]:
    M_opt, ff_opt = opt_mach_for_ci(CI)
    print(f"{CI:14} {M_opt:12.3f} {ff_opt:.1f}")

# Plot: fuel flow and SR vs Mach
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(M_range, ff_kghr, color="#2563EB", lw=2)
ax1.set_xlabel("Mach");  ax1.set_ylabel("Fuel Flow (kg/hr)")
ax1.set_title("Fuel Flow vs Mach");  ax1.grid(True, alpha=0.3)
ax2.plot(M_range, SR,      color="#16A34A", lw=2)
ax2.set_xlabel("Mach");  ax2.set_ylabel("Specific Range (km/kg)")
ax2.set_title("Specific Range vs Mach");  ax2.grid(True, alpha=0.3)
plt.tight_layout();  plt.show()

import numpy as np

def breguet_range(V_ms, LD, SFC_kgNs, W_initial, W_final):
    """
    Breguet range equation.
    V_ms      : cruise TAS, m/s
    LD        : cruise lift-to-drag ratio
    SFC_kgNs  : specific fuel consumption, kg/(N·s)
    W_initial : initial cruise weight, N
    W_final   : final cruise weight (= W_initial - W_fuel), N
    Returns   : range in km
    """
    return (V_ms / (SFC_kgNs * 9.80665)) * LD * np.log(W_initial / W_final) / 1000

# A320neo-like parameters
MTOW    = 79000    # kg — max takeoff weight
OEW     = 42600    # kg — operating empty weight
payload = 15000    # kg — passengers + bags (~150 pax)
fuel    = 18000    # kg — usable fuel

h_m     = 35000 * 0.3048
_, _, _, a = isa(h_m)
V_cruise= 0.78 * a         # TAS at M0.78, FL350
LD      = 17.5              # typical A320neo cruise L/D
SFC     = 1.65e-5           # kg/(N·s) — LEAP-1A26 typical

W_init  = (OEW + payload + fuel) * 9.80665
W_final = (OEW + payload) * 9.80665

R = breguet_range(V_cruise, LD, SFC, W_init, W_final)
print(f"Cruise Mach:     M0.78")
print(f"TAS:             {V_cruise/0.514444:.1f} kts")
print(f"L/D:             {LD}")
print(f"Fuel burned:     {fuel:,} kg")
print(f"Range (Breguet): {R:.0f} km  ({R/1.852:.0f} nm)")

# Sensitivity: range vs payload (range-payload curve)
payloads  = np.linspace(0, 20000, 100)
ranges    = []
for pl in payloads:
    fuel_avail = MTOW - OEW - pl        # fuel limited by MTOW
    fuel_avail = np.clip(fuel_avail, 0, 18000)
    Wi = (OEW + pl + fuel_avail) * 9.80665
    Wf = (OEW + pl) * 9.80665
    ranges.append(breguet_range(V_cruise, LD, SFC, Wi, Wf))

import matplotlib.pyplot as plt
plt.figure(figsize=(9, 5))
plt.plot(payloads/1000, ranges, color="#2563EB", lw=2)
plt.xlabel("Payload (tonnes)");  plt.ylabel("Range (km)")
plt.title("Range–Payload Curve (A320neo-like)")
plt.grid(True, alpha=0.3);  plt.tight_layout();  plt.show()

import numpy as np
from scipy.integrate import solve_ivp

# Aircraft (A320-like at typical climb weight)
W_N   = 70000 * 9.80665   # N
S     = 122.6               # m²
CD0   = 0.0280              # climb config (no flaps)
k     = 0.0375
SFC   = 1.80e-5             # kg/(N·s) climb SFC (higher than cruise)

# Thrust model: T = T_sl * (rho/rho0)^0.75  (simplified)
T_SL  = 120000              # N  — total SL static thrust (2 × CFM56)

def climb_ode(t, y):
    """y = [altitude m, fuel burned kg]"""
    h, mf = y
    T_ISA, P, rho, a = isa(h)
    T    = T_SL * (rho / 1.225) ** 0.75      # thrust at altitude
    M    = 0.76                                 # constant Mach climb above 10000ft
    V    = M * a
    q    = 0.5 * rho * V**2
    CL   = W_N / (q * S)
    CD   = CD0 + k * CL**2
    D    = q * S * CD
    ROC  = V * (T - D) / W_N               # m/s rate of climb
    ff   = T * SFC                             # kg/s fuel flow
    return [ROC, ff]

# Integrate from 3000ft to FL350
h_init = 3000 * 0.3048
h_end  = 35000 * 0.3048

# Stop event: when altitude reaches FL350
def reached_cruise(t, y): return y[0] - h_end
reached_cruise.terminal = True

sol = solve_ivp(climb_ode, [0, 3600], [h_init, 0],
                events=reached_cruise, max_step=10,
                dense_output=True)

t_climb    = sol.t_events[0][0]
fuel_climb = sol.y_events[0][0][1]

print(f"Time to climb FL030→FL350: {t_climb/60:.1f} min")
print(f"Fuel to climb:             {fuel_climb:.0f} kg")
print(f"Average ROC:               {(h_end-h_init)/t_climb*196.85:.0f} fpm")

import numpy as np

def opt_altitude(W_kg, M=0.78):
    """
    Find altitude that maximises specific range for given weight and Mach.
    Uses simplified drag model.
    """
    altitudes = np.linspace(25000, 43000, 200) * 0.3048
    best_SR, best_alt = 0, altitudes[0]
    for h in altitudes:
        T, P, rho, a = isa(h)
        V  = M * a
        q  = 0.5 * rho * V**2
        CL = W_kg * 9.80665 / (q * 122.6)
        CD = 0.0240 + 0.0375 * CL**2
        D  = q * 122.6 * CD
        ff = D * V / (0.30 * 43.2e6)
        SR = V / ff / 1000
        if SR > best_SR:
            best_SR, best_alt = SR, h
    return best_alt / 0.3048, best_SR

# Track optimum altitude as fuel burns off during flight
weights = np.arange(75000, 52000, -1000)   # kg, fuel burning off
print(f"{'Weight (kg)':14} {'Opt Alt (ft)':15} {'Spec. Range (km/kg)'}")
print("-"*46)
for W in weights[::5]:           # print every 5th step
    alt, sr = opt_altitude(W)
    print(f"{W:14,} {alt:15.0f} {sr:.4f}")

import numpy as np
import matplotlib.pyplot as plt

# Aircraft parameters
W     = 70000    # kg  — design weight
S     = 122.6    # m²
CLmax = 1.52     # clean configuration
CLmin = -0.80    # negative lift limit

# CS-25 load factor limits
n_lim_pos = 2.5   # positive limit (CS25.337)
n_lim_neg = -1.0  # negative limit
n_ult_pos = n_lim_pos * 1.5   # 1.5 × limit = ultimate

# Reference conditions (sea level ISA)
rho   = 1.225
g     = 9.80665
W_N   = W * g

# Stall speed at 1g (EAS)
Vs    = np.sqrt(2 * W_N / (rho * S * CLmax))   # m/s EAS
Vs_kts= Vs / 0.514444

# Manoeuvre speed Va = Vs * sqrt(n_lim)
Va    = Vs * np.sqrt(n_lim_pos)
Vc    = 175        # m/s EAS — design cruise speed
Vd    = 210        # m/s EAS — design dive speed

# Stall boundary curves
V_range = np.linspace(0, Vd, 500)
n_pos   = (rho * V_range**2 * S * CLmax) / (2 * W_N)   # positive stall
n_neg   = (rho * V_range**2 * S * CLmin) / (2 * W_N)   # negative stall

fig, ax = plt.subplots(figsize=(10, 6))

# Positive stall boundary (up to n_lim)
mask_pos = n_pos <= n_lim_pos
ax.plot(V_range[mask_pos] / 0.514444, n_pos[mask_pos],
        color="#2563EB", lw=2, label="Stall boundary")

# Structural limits
ax.axhline(n_lim_pos, color="#DC2626", lw=1.5, ls="--", label=f"Limit load n={n_lim_pos}g")
ax.axhline(n_ult_pos, color="#DC2626", lw=1,   ls=":",  label=f"Ultimate load n={n_ult_pos}g")
ax.axhline(n_lim_neg, color="#2563EB", lw=1.5, ls="--", label=f"Negative limit n={n_lim_neg}g")

# Vertical speed limits
ax.axvline(Va  / 0.514444, color="gray", lw=1, ls=":", label=f"Va = {Va/0.514444:.0f} kts")
ax.axvline(Vc  / 0.514444, color="green", lw=1, ls=":", label=f"Vc = {Vc/0.514444:.0f} kts")
ax.axvline(Vd  / 0.514444, color="red",   lw=1, ls=":", label=f"Vd = {Vd/0.514444:.0f} kts")

ax.fill_between(V_range[mask_pos] / 0.514444,
                 n_pos[mask_pos], 0, alpha=0.07, color="#2563EB")
ax.set_xlabel("EAS (knots)");  ax.set_ylabel("Load Factor n (g)")
ax.set_title("V-n Diagram — CS-25 Transport")
ax.legend(fontsize=9);  ax.grid(True, alpha=0.3)
ax.set_xlim(0, Vd/0.514444 + 20);  ax.set_ylim(-1.5, 4.5)
plt.tight_layout();  plt.show()
print(f"1g stall speed: {Vs_kts:.1f} kts EAS")
print(f"Va (manoeuvre): {Va/0.514444:.1f} kts EAS")

import numpy as np

# S-N curve for Al 2024-T3 (log-log linear fit from AV3 regression)
# log10(N) = a - b * log10(sigma)
a_sn = 15.2
b_sn = 4.8

def cycles_to_failure(sigma_MPa):
    """N from S-N curve"""
    return 10 ** (a_sn - b_sn * np.log10(sigma_MPa))

# Mission spectrum: typical short-haul aircraft
# Each row: [stress amplitude MPa, occurrences per flight]
spectrum = np.array([
    [250,  1   ],   # rotation / landing
    [180,  3   ],   # turbulence moderate
    [120,  8   ],   # light turbulence
    [80,   25  ],   # manoeuvres
    [50,   100 ],   # pressurisation cycles within flight
    [30,   500 ],   # minor vibration
])

sigma = spectrum[:, 0]
n_occ = spectrum[:, 1]
N_f   = cycles_to_failure(sigma)
D_per_flight = np.sum(n_occ / N_f)   # Miner's sum per flight
life_flights  = 1 / D_per_flight       # flights until D = 1.0

print(f"{'Stress (MPa)':15} {'n/flight':12} {'N_f':15} {'n/N (damage)'}")
print("-"*55)
for i in range(len(sigma)):
    print(f"{sigma[i]:15.0f} {n_occ[i]:12.0f} {N_f[i]:15.2e} {n_occ[i]/N_f[i]:.2e}")
print(f"\nDamage per flight: {D_per_flight:.2e}")
print(f"Predicted life:    {life_flights:,.0f} flights")
print(f"                   {life_flights/365:.0f} years at 1 flight/day")

# Current fleet status
flights_flown   = 12500
damage_accrued  = flights_flown * D_per_flight
remaining_life  = (1 - damage_accrued) / D_per_flight
print(f"\nFlights flown:     {flights_flown:,}")
print(f"Damage accrued:    {damage_accrued:.3f}  ({damage_accrued*100:.1f}%)")
print(f"Remaining life:    {remaining_life:,.0f} flights")

import numpy as np
import matplotlib.pyplot as plt

# Flight test data: (CL, CD) pairs from level-flight test points
# Each row is one test condition (different weight/altitude/speed)
CL_data = np.array([0.20, 0.30, 0.40, 0.50, 0.60,
                     0.70, 0.80, 0.90, 1.00, 1.10])
CD_data = np.array([0.0258, 0.0268, 0.0284, 0.0311, 0.0347,
                     0.0393, 0.0454, 0.0529, 0.0625, 0.0741])

# Fit CD = CD0 + k * CL^2  →  linear regression on (CL², CD)
CL2      = CL_data**2
coeffs   = np.polyfit(CL2, CD_data, 1)
k_fit    = coeffs[0]
CD0_fit  = coeffs[1]

# Oswald efficiency factor
AR       = 9.4         # A320 aspect ratio
e_oswald = 1 / (np.pi * AR * k_fit)

print(f"Drag polar fit:")
print(f"  CD0   = {CD0_fit:.5f}")
print(f"  k     = {k_fit:.5f}")
print(f"  e     = {e_oswald:.3f}  (Oswald efficiency)")

# Maximum L/D and the CL at which it occurs
CL_maxLD = np.sqrt(CD0_fit / k_fit)
LD_max   = CL_maxLD / (2 * np.sqrt(CD0_fit * k_fit))
print(f"  CL at max L/D = {CL_maxLD:.3f}")
print(f"  Max L/D       = {LD_max:.1f}")

# Plot the polar
CL_curve = np.linspace(0, 1.3, 200)
CD_curve = CD0_fit + k_fit * CL_curve**2

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# CL vs CD
ax1.scatter(CD_data, CL_data, color="steelblue", s=50, zorder=5, label="Test data")
ax1.plot(CD_curve, CL_curve, "r-", lw=2, label=f"Fit: CD0={CD0_fit:.4f}, k={k_fit:.4f}")
ax1.axvline(2*np.sqrt(CD0_fit*k_fit), color="green", ls="--", lw=1, label=f"Min CD (max L/D)")
ax1.set_xlabel("CD");  ax1.set_ylabel("CL")
ax1.set_title("Drag Polar");  ax1.legend(fontsize=8);  ax1.grid(True, alpha=0.3)

# L/D vs CL
LD_curve = CL_curve / CD_curve
ax2.plot(CL_curve, LD_curve, color="#16A34A", lw=2)
ax2.axvline(CL_maxLD, color="red", ls="--", lw=1, label=f"Best CL={CL_maxLD:.3f}, L/D={LD_max:.1f}")
ax2.set_xlabel("CL");  ax2.set_ylabel("L/D")
ax2.set_title("Lift-to-Drag Ratio");  ax2.legend(fontsize=9);  ax2.grid(True, alpha=0.3)
plt.tight_layout();  plt.show()

import pandas as pd
import numpy  as np
import matplotlib.pyplot as plt

# QAR data is usually a CSV or binary file exported from the ACMS
# Typical columns: time, altitude, airspeed, heading, N1%, temp,
#                  fuel_flow, load_factor, pitch, roll, etc.
df = pd.read_csv("flight_WB123_20260322.csv")

# Basic inspection
print(f"Flight duration: {len(df)/4:.0f} seconds ({len(df)/4/60:.1f} min)")  # 4 Hz data
print(f"Parameters: {list(df.columns)}")
print(df.describe().round(2))

# Check for invalid/missing data
print("\nMissing values:")
print(df.isnull().sum())

# Replace clearly invalid sensor spikes (e.g. altitude = 0 mid-flight)
df["altitude_ft"] = df["altitude_ft"].replace(0, np.nan).interpolate()

import pandas as pd
import numpy  as np

# Define exceedance thresholds (simplified A320 limits)
limits = {
    "load_factor_nz":  (-0.5,   2.0 ),   # g  — smooth air ops limits
    "pitch_deg":       (-5.0,   25.0),   # deg
    "bank_angle_deg":  (-30.0,  30.0),   # deg
    "N1_pct":          (0,      104.0),  # %
    "airspeed_kts":    (0,      350.0),  # kts below FL100
    "vsi_fpm":         (-2000, 2000 ),   # fpm vertical speed
}

exceedances = []

for param, (lo, hi) in limits.items():
    if param not in df.columns:
        continue
    col   = df[param]
    exceed= (col < lo) | (col > hi)

    if exceed.any():
        # Find contiguous exceedance windows
        changes = exceed.astype(int).diff().fillna(0)
        starts  = df.index[changes == 1]
        ends    = df.index[changes == -1]
        for s, e in zip(starts, ends):
            window   = col.loc[s:e]
            peak_val = window.max() if window.max() > hi else window.min()
            duration = (e - s) / 4              # seconds (4Hz data)
            exceedances.append({
                "parameter": param,
                "limit":     hi if peak_val > hi else lo,
                "peak_value": peak_val.round(2),
                "duration_s": duration,
                "time_s":     s / 4
            })

result = pd.DataFrame(exceedances)
if result.empty():
    print("No exceedances detected. Clean flight.")
else:
    print(f"{len(result)} exceedances detected:")
    print(result.to_string(index=False))

import pandas as pd
import numpy  as np
import matplotlib.pyplot as plt

# Load month of flight records (from OPS system export)
flights = pd.read_csv("fleet_march2026.csv")
# Expected columns: flight_no, reg, route, dist_nm, fuel_plan_kg,
#                   fuel_actual_kg, pax, cargo_kg, wind_comp_kts

# Fuel delta: actual minus planned
flights["fuel_delta"]     = flights["fuel_actual_kg"] - flights["fuel_plan_kg"]
flights["fuel_delta_pct"] = flights["fuel_delta"] / flights["fuel_plan_kg"] * 100
flights["fuel_per_nm"]    = flights["fuel_actual_kg"] / flights["dist_nm"]

# Fleet summary
print("=== Fleet Fuel Performance — March 2026 ===")
print(f"Flights analysed:     {len(flights)}")
print(f"Mean fuel delta:      {flights.fuel_delta.mean():.0f} kg ({flights.fuel_delta_pct.mean():.1f}%)")
print(f"Std deviation:        {flights.fuel_delta.std():.0f} kg")
print(f"Worst flight:         +{flights.fuel_delta.max():.0f} kg")
print(f"Best flight:          {flights.fuel_delta.min():.0f} kg")

# Per-tail analysis (identify aircraft with chronic over-burn)
per_reg = flights.groupby("reg").agg(
    flights       = ("flight_no", "count"),
    mean_delta_kg = ("fuel_delta", "mean"),
    total_excess  = ("fuel_delta", "sum")
).round(1).sort_values("mean_delta_kg", ascending=False)
print("\nPer-aircraft fuel performance:")
print(per_reg.head(10))

# Correlation: wind vs fuel delta
corr = flights["wind_comp_kts"].corr(flights["fuel_delta"])
print(f"\nWind vs fuel delta correlation: {corr:.3f}")
# Strong negative correlation expected: headwind → more fuel

# Plot distribution of fuel delta
plt.figure(figsize=(9, 4))
plt.hist(flights["fuel_delta"], bins=40, color="steelblue", edgecolor="white")
plt.axvline(0, color="red", lw=1.5, label="Plan")
plt.axvline(flights.fuel_delta.mean(), color="orange", lw=1.5, ls="--", label=f"Mean: {flights.fuel_delta.mean():.0f} kg")
plt.xlabel("Fuel Delta (kg)");  plt.ylabel("Flights")
plt.title("Actual vs Planned Fuel — March 2026")
plt.legend();  plt.tight_layout();  plt.show()