#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Richards S-Curve Fit for 1984-2025 Data Fits and visualizes the Richards growth model to rocket launch data starting from 1984. """ import pandas as pd import numpy as np from scipy.optimize import curve_fit import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import warnings warnings.filterwarnings('ignore') plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] plt.rcParams['axes.unicode_minus'] = False def richards(t, K, r, t0, v): """ Richards curve (generalized logistic model) N(t) = K / (1 + exp(-r*(t - t0)))^(1/v) """ exp_term = np.exp(-r * (t - t0)) exp_term = np.clip(exp_term, 1e-10, 1e10) return K / np.power(1 + exp_term, 1/v) def load_data(filepath="rocket_launch_counts.csv"): """Load rocket launch data""" df = pd.read_csv(filepath) df = df.rename(columns={"YDate": "year", "Total": "launches"}) df["year"] = pd.to_numeric(df["year"], errors="coerce") df["launches"] = pd.to_numeric(df["launches"], errors="coerce") df = df.dropna(subset=["year", "launches"]) df = df[(df["year"] >= 1957) & (df["year"] <= 2025)] df = df.astype({"year": int, "launches": int}) df = df.sort_values("year").reset_index(drop=True) return df def fit_richards(years, launches, base_year=1984): """Fit Richards model to data""" t = (years - base_year).astype(float) # Initial parameters and bounds (unconstrained K) p0 = [5000.0, 0.1, 40.0, 2.0] bounds = ([500, 0.005, 10, 0.2], [100000, 1.0, 150, 10.0]) popt, pcov = curve_fit(richards, t, launches, p0=p0, bounds=bounds, maxfev=100000) perr = np.sqrt(np.diag(pcov)) # Calculate R² y_pred = richards(t, *popt) ss_res = np.sum((launches - y_pred) ** 2) ss_tot = np.sum((launches - np.mean(launches)) ** 2) r_squared = 1 - (ss_res / ss_tot) return popt, perr, r_squared def main(): print("=" * 60) print("Richards S-Curve Fit (1984-2025)") print("=" * 60) # Load data df = load_data() # Filter 1984-2025 start_year = 1984 data = df[(df["year"] >= start_year) & (df["year"] <= 2025)].copy() years = data["year"].values launches = data["launches"].values print(f"Data range: {start_year} - 2025") print(f"Data points: {len(data)}") # Fit model popt, perr, r_squared = fit_richards(years, launches, base_year=start_year) K, r, t0, v = popt print(f"\nFitted Parameters:") print(f" K (carrying capacity) = {K:.0f} launches/year") print(f" r (growth rate) = {r:.4f}") print(f" t0 (inflection point) = {start_year + t0:.1f}") print(f" v (shape parameter) = {v:.3f}") print(f" R² = {r_squared:.4f}") # Physical limit physical_max = 3650 print(f"\nPhysical limit: {physical_max} (10 sites × 365 days)") print(f"K / Physical limit = {K/physical_max:.2f}x") # ========== Create Visualization ========== fig, ax = plt.subplots(figsize=(12, 7)) # Historical data points ax.scatter(years, launches, color='#2C3E50', s=80, alpha=0.8, label='Historical Data (1984-2025)', zorder=3, edgecolor='white', linewidth=0.5) # Generate smooth S-curve years_smooth = np.linspace(start_year, 2100, 500) t_smooth = years_smooth - start_year pred_smooth = richards(t_smooth, *popt) # S-curve prediction ax.plot(years_smooth, pred_smooth, color='#27AE60', lw=3, label=f'Richards Model (K={K:.0f}, R²={r_squared:.3f})', zorder=2) # Confidence band (approximate using parameter errors) K_low = max(500, K - 2*perr[0]) K_high = K + 2*perr[0] pred_low = richards(t_smooth, K_low, r, t0, v) pred_high = richards(t_smooth, K_high, r, t0, v) ax.fill_between(years_smooth, pred_low, pred_high, color='#27AE60', alpha=0.15, label='95% Confidence Band') # Physical limit line ax.axhline(physical_max, color='#E74C3C', ls='--', lw=2.5, label=f'Physical Limit: {physical_max} (1/site/day)') # Mark inflection point inflection_year = start_year + t0 inflection_value = K / (v + 1) ** (1/v) ax.scatter([inflection_year], [inflection_value], color='#F39C12', s=150, marker='*', zorder=5, label=f'Inflection Point ({inflection_year:.0f})') # Mark key years ax.axvline(2025, color='gray', ls=':', lw=1.5, alpha=0.7) ax.axvline(2050, color='#3498DB', ls=':', lw=2, alpha=0.8) ax.text(2026, K*0.95, '2025\n(Now)', fontsize=10, color='gray') ax.text(2051, K*0.85, '2050\n(Target)', fontsize=10, color='#3498DB') # Prediction for 2050 t_2050 = 2050 - start_year pred_2050 = richards(t_2050, *popt) ax.scatter([2050], [pred_2050], color='#3498DB', s=100, marker='D', zorder=4) ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), xytext=(2055, pred_2050-300), fontsize=10, color='#3498DB', fontweight='bold') # Formatting ax.set_xlabel('Year', fontsize=12) ax.set_ylabel('Annual Launches', fontsize=12) ax.set_title('Richards S-Curve Model Fit (1984-2025 Data)\nRocket Launch Capacity Projection', fontsize=14, fontweight='bold') ax.legend(loc='upper left', fontsize=10) ax.grid(True, alpha=0.3) ax.set_xlim(1982, 2100) ax.set_ylim(0, K * 1.15) # Add text box with model info textstr = f'''Model: N(t) = K / (1 + e^(-r(t-t₀)))^(1/v) Data Window: {start_year}-2025 ({len(data)} points) K = {K:.0f} launches/year r = {r:.4f} t₀ = {start_year + t0:.1f} v = {v:.3f} R² = {r_squared:.4f} Note: Physical limit (3,650) shown as dashed red line''' props = dict(boxstyle='round', facecolor='wheat', alpha=0.8) ax.text(0.98, 0.35, textstr, transform=ax.transAxes, fontsize=9, verticalalignment='top', horizontalalignment='right', bbox=props, family='monospace') plt.tight_layout() plt.savefig('richards_curve_1984.png', dpi=150, bbox_inches='tight') plt.close() print("\nPlot saved: richards_curve_1984.png") print("=" * 60) if __name__ == "__main__": main()